Hello, everyone. Recently, I have been solving the compressible NS equations on the unit disk, with the model as follows:
where,
J =E + u \times B.
and
The numerical scheme is：
When solving this equation, the pressures p, electric fields E and magnetic fields B all involved were substituted using the exact solutions. If the functions ρ and u are discretized using the P1 element, the L^2 convergence order obtained in both cases is 2. Now my main goal is to test the calculation of higher-order elements in curved domain.
However, when I discretize using the P2 element, the L2 error convergence order of rho can only reach 2nd order, while the L2 error convergence order of u can reach 3rd order.
Here is my code：
import gmsh
from dolfinx.io import gmshio
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.fem import Function
import time
from ufl import dx, grad, dot,inner, div,nabla_grad,curl,split, lhs, rhs
from basix.ufl import element,mixed_element
from numpy import sin, cos, pi, exp
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
# 精确解
mu = 1.0
muc = 2.0
Rm=1.0
gamma=5/3
kappa=1.0
def ExactRho(x, t): # 用恒正的 ρ 能解决问题——所有内部节点的 ρ 都远离零，矩阵条件数好
r2 = x[0]**2 + x[1]**2
return (exp(-t) + 1.0) * (1 - r2)
def ExactU(x,t):
r2 = x[0]**2 + x[1]**2
factor = t * (1 - r2)
u1 = factor * sin(pi*x[0]) * cos(pi*x[1])
u2 = factor * cos(pi*x[0]) * sin(pi*x[1])
return [u1, u2]
def Exactp(x,t):
r2 = x[0]**2 + x[1]**2
return (exp(-t) + 1.0) * (1 - r2) * sin(pi*x[0])
def ExactE(x,t):
r2 = x[0]**2 + x[1]**2
fr = 1 - r2
return pi * fr**2 * sin(pi*x[0]) * sin(pi*x[1]) * sin(pi*t)
def ExactB(x, t):
xv, yv = x[0], x[1]
r2 = xv**2 + yv**2
fr = 1 - r2
phi_f = cos(pi*t) + 2.0
sinx = sin(pi*xv); siny = sin(pi*yv)
cosx = cos(pi*xv); cosy = cos(pi*yv)
B1 = fr * phi_f * sinx * (pi*fr*cosy - 4*yv*siny)
B2 = -fr * phi_f * siny * (pi*fr*cosx - 4*xv*sinx)
return [B1, B2]
def rhs_expr0(x, t):
xv, yv = x[0], x[1]
r2 = xv**2 + yv**2
fr = 1 - r2
sinx = sin(pi*xv); cosx = cos(pi*xv)
siny = sin(pi*yv); cosy = cos(pi*yv)
exp_t = exp(-t)
f1 = 1.0 + exp_t
h = (- 4*t*xv*f1*fr*sinx*cosy
- 4*t*yv*f1*fr*siny*cosx
+ 2*pi*t*f1*fr*fr*cosx*cosy
- fr*exp_t)
return h
def rhs_expr1(x, t):
xv, yv = x[0], x[1]
r2 = xv**2 + yv**2
fr = 1 - r2
sinx = sin(pi*xv)
cosx = cos(pi*xv)
siny = sin(pi*yv)
cosy = cos(pi*yv)
sint = sin(pi*t)
phi_f = cos(pi*t) + 2.0
# ρ
rho = (exp(-t) + 1.0) * (1 - r2)
# u
u1 = t * fr * sinx * cosy
u2 = t * fr * cosx * siny
# ∂u/∂t
du1_dt = fr * sinx * cosy
du2_dt = fr * cosx * siny
# ∂u1/∂x, ∂u1/∂y
du1_dx = t * (-2*xv*sinx*cosy + fr*pi*cosx*cosy)
du1_dy = t * (-2*yv*sinx*cosy - fr*pi*sinx*siny)
# ∂u2/∂x, ∂u2/∂y
du2_dx = t * (-2*xv*cosx*siny - fr*pi*sinx*siny)
du2_dy = t * (-2*yv*cosx*siny + fr*pi*cosx*cosy)
# u·∇u
u_grad_u1 = u1*du1_dx + u2*du1_dy
u_grad_u2 = u1*du2_dx + u2*du2_dy
# Δu
# ∂²u1/∂x²
d2u1_dx2 = t * (-2*sinx*cosy - 4*xv*pi*cosx*cosy - fr*pi**2*sinx*cosy)
# ∂²u1/∂y²
d2u1_dy2 = t * (-2*sinx*cosy + 4*yv*pi*sinx*siny - fr*pi**2*sinx*cosy)
lap_u1 = d2u1_dx2 + d2u1_dy2
# ∂²u2/∂x²
d2u2_dx2 = t * (-2*cosx*siny + 4*xv*pi*sinx*siny - fr*pi**2*cosx*siny)
# ∂²u2/∂y²
d2u2_dy2 = t * (-2*cosx*siny - 4*yv*pi*cosx*cosy - fr*pi**2*cosx*siny)
lap_u2 = d2u2_dx2 + d2u2_dy2
# \nabla p
dp_dx = (exp(-t) + 1.0)* (pi*fr*cosx - 2*xv*sinx)
dp_dy = (exp(-t) + 1.0) * (-2*yv*sinx)
# E和B
E_val = pi * fr*fr * sinx * siny * sint
Px = pi*fr*cosx - 4*xv*sinx
Py = pi*fr*cosy - 4*yv*siny
B1 = fr * phi_f * sinx * Py
B2 = -fr * phi_f * siny * Px
# E×B = [-E*B2, E*B1]
E_cross_B1 = -E_val * B2
E_cross_B2 = E_val * B1
# u×B×B = [-(u×B)*B2, (u×B)*B1]
u_cross_B = u1*B2 - u2*B1
u_B_B1 = -u_cross_B * B2
u_B_B2 = u_cross_B * B1
# ∇(∇·u)
d_div_u_dx = t * (-2*sinx*cosy - 6*xv*pi*cosx*cosy + 2*yv*pi*sinx*siny - 2*fr*pi*pi*sinx*cosy)
d_div_u_dy = t * (2*xv*pi*sinx*siny - 2*cosx*siny - 6*yv*pi*cosx*cosy - 2*fr*pi*pi*cosx*siny)
# 组装 f = ρ*(∂u/∂t + u·∇u) + \nabla p-J\times B- Δu
f1 = rho*du1_dt + rho*u_grad_u1 + dp_dx - E_cross_B1-u_B_B1-(2*muc-mu)*d_div_u_dx- mu*lap_u1
f2 = rho*du2_dt + rho*u_grad_u2+ dp_dy - E_cross_B2-u_B_B2-(2*muc-mu)*d_div_u_dy- mu*lap_u2
return [f1, f2]
def cross_2d(a, b):
return a[0]*b[1] - a[1]*b[0]
def cross1(a, b):
return ufl.as_vector([-a*b[1], a*b[0]])
def cross1_rev(b, a):
return ufl.as_vector([a * b[1], -a * b[0]])
def clamped_boundary(x):
return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)
# 固定终止时间
T = 1
# 结果记录
hs = []
errorsrho_L2 = []
errorsrho_H1 = []
errorsU_L2 = []
errorsU_H1 = []
start_time = time.time()
for i in range(3):
n = 2 ** (i + 2)
h = 1.0 / n
num_steps = 2 ** (3 * i + 6)
dt = T / num_steps
t = 0.0
# 生成网格
gmsh.initialize()
membrane = gmsh.model.occ.addDisk(0,0,0,1,1)
gmsh.model.occ.synchronize()
gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)
gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.model.mesh.generate(gdim)
domain, _, _ = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)
gmsh.finalize()
x = ufl.SpatialCoordinate(domain)
V = fem.functionspace(domain, ("Lagrange", 2)) # 线性 Lagrange 元
VE_element = element("Lagrange", domain.basix_cell(), 2)
VB_element = element("RT", domain.basix_cell(), 2)
VU_element = element("Lagrange", domain.basix_cell(), 2, shape=(2, ))
V_element = element("Lagrange", domain.basix_cell(), 2)
Vv = fem.functionspace(domain, VU_element)
VE = fem.functionspace(domain, VE_element)
VB = fem.functionspace(domain, VB_element)
VM = fem.functionspace(domain, mixed_element([VU_element, V_element])) # 混合有限元空间
V0 = VM.sub(0)
V1 = VM.sub(1)
Q, _ = V0.collapse() # u 平铺成了标量空间
S, _ = V1.collapse() # rho
## 从这里开始计算
u_n = fem.Function(Q)
u_n.interpolate(lambda x:ExactU(x,t))
p_n = fem.Function(V)
p_n.interpolate(lambda x: Exactp(x,t))
rho_n = fem.Function(S)
rho_n.interpolate(lambda x: ExactRho(x, t))
# 边界条件
rho_D = fem.Function(S)
rho_D.interpolate(lambda x: ExactRho(x, t))
u_D = fem.Function(Q)
u_D.interpolate(lambda x: ExactU(x, t))
B_n = fem.Function(VB)
B_n.interpolate(lambda x: ExactB(x, t))
E_n = fem.Function(VE)
# 右端函数定义
f_fun0 = fem.Function(S)
f_fun1 = fem.Function(Q)
f_fun0.interpolate(lambda x: rhs_expr0(x, t))
f_fun1.interpolate(lambda x: rhs_expr1(x, t))
# 施加 Dirichlet 边界条件（3D 边界）
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
# 施加 Dirichlet 边界条件
facets = mesh.locate_entities_boundary(
domain,
dim=(domain.topology.dim - 1),
marker=clamped_boundary)
boundary_dofs0 = fem.locate_dofs_topological((V0, Q), entity_dim=1, entities=facets)
boundary_dofs3 = fem.locate_dofs_topological((V1, S), entity_dim=1, entities=facets)
bc0 = fem.dirichletbc(value=u_D, dofs=boundary_dofs0, V=V0)
bc3 = fem.dirichletbc(value=rho_D, dofs=boundary_dofs3, V=V1)
bc=[bc0,bc3]
# 在时间循环前创建 wh_n 并设置初始条件
wh_n = Function(VM)
# 设置初始条件到 wh_n
u_init = fem.Function(Q)
u_init.interpolate(lambda x: ExactU(x, 0.0)) # t=0的初始条件
rho_init = fem.Function(S)
rho_init.interpolate(lambda x: ExactRho(x, 0.0))
wh_n.sub(0).interpolate(u_init)
wh_n.sub(1).interpolate(rho_init)
# u E B 的变分形式
wh = Function(VM)
u_h = fem.Function(Q)
rho_h = fem.Function(S)
u, rho = split(wh)
v, w = ufl.TestFunctions(VM)
F1 = rho * w * dx - rho_n * w * dx - dt * dot(rho * u, grad(w)) * dx - dt *f_fun0* w * dx \
+ rho_n* dot(u, v) *dx - rho_n* dot(u_n, v) *dx\
+ dt*(1/2)*inner(rho_n*grad(dot(u,u)),v)*dx - dt*inner(rho_n*cross1_rev(u,curl(u_n)), v)*dx \
+ dt* mu * inner(nabla_grad(u),nabla_grad(v))*dx + dt*(2*muc-mu)*inner(div(u),div(v))*dx- dt* inner(f_fun1,v)*dx\
- dt* dot(p_n, div(v))*dx - dt*inner(cross1(cross_2d(u,B_n),B_n), v) * dx-dt*inner(cross1(E_n, B_n), v) * dx
for n_step in range(num_steps):
t += dt
f_fun0.interpolate(lambda x: rhs_expr0(x, t))
f_fun1.interpolate(lambda x: rhs_expr1(x, t))
E_n.interpolate(lambda x:ExactE(x, t))
wh.x.array[:] = wh_n.x.array[:]
problem = NonlinearProblem(F1, wh, bcs=bc)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-10
solver.rtol = 1e-10
solver.max_it = 100
solver.krylov_options = {
"ksp_type": "preonly",
"pc_type": "lu",
}
n_iter, converged = solver.solve(wh)
u_h = wh.sub(0).collapse()
rho_h = wh.sub(1).collapse()
if not converged:
print(f" WARNING: N={n}, step={num_steps+1}, "
f"t={t:.4f}: not converged ({n_iter} iter)")
wh_n.x.array[:] = wh.x.array[:]
rho_n.x.array[:] = rho_h.x.array
u_n.x.array[:] = u_h.x.array
p_n.interpolate(lambda x:Exactp(x, t))
B_n.interpolate(lambda x:ExactB(x, t))
# 计算误差
V_ex = fem.functionspace(domain, ("Lagrange", 3))
rho_ex = fem.Function(V_ex)
rho_ex.interpolate(lambda x: ExactRho(x,T))
# E 的误差
E_ex = fem.Function(V_ex)
E_ex.interpolate(lambda x: ExactE(x,T))
p_ex = fem.Function(V_ex)
p_ex.interpolate(lambda x: Exactp(x,T))
Vv_ex = fem.functionspace(domain, ("Lagrange", 3, (2, )))
# u 的误差
u_ex = fem.Function(Vv_ex)
u_ex.interpolate(lambda x: ExactU(x,T))
# B 的误差
B_ex = fem.Function(Vv_ex)
B_ex.interpolate(lambda x: ExactB(x,T))
# L2 误差
errorrho_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((rho_h - rho_ex)**2 * ufl.dx)), op=MPI.SUM))
errorrho_H1 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form(dot(grad(rho_h - rho_ex),grad(rho_h - rho_ex)) * ufl.dx)), op=MPI.SUM))
errorU_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((u_h - u_ex)**2 * ufl.dx)), op=MPI.SUM))
errorU_H1 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form(inner(grad(u_h - u_ex),grad(u_h - u_ex)) * ufl.dx)), op=MPI.SUM))
if MPI.COMM_WORLD.rank == 0:
hs.append(h)
errorsrho_L2.append(errorrho_L2)
errorsrho_H1.append(errorrho_H1)
errorsU_L2.append(errorU_L2)
errorsU_H1.append(errorU_H1)
print(f"[i={i}] n={n}, num_steps={num_steps}, dt={dt:.3e}, h={h:.3e}")
print(f" → L2 errorrho = {errorrho_L2:.3e}, H1 errorrho = {errorrho_H1:.3e}")
print(f" → L2 errorU = {errorU_L2:.3e}, H1 errorU = {errorU_H1:.3e}") #
end_time = time.time() # 计时结束
elapsed = end_time - start_time
if MPI.COMM_WORLD.rank == 0:
print("\nEstimated convergence rates (L2):")
for j in range(1, len(hs)):
e1 = np.log(errorsrho_L2[j-1]/errorsrho_L2[j]) / np.log(hs[j-1]/hs[j])
e2 = np.log(errorsrho_H1[j-1]/errorsrho_H1[j]) / np.log(hs[j-1]/hs[j])
r1 = np.log(errorsU_L2[j-1]/errorsU_L2[j]) / np.log(hs[j-1]/hs[j])
r2 = np.log(errorsU_H1[j-1]/errorsU_H1[j]) / np.log(hs[j-1]/hs[j])
print(f"h={hs[j-1]:.3e} → h={hs[j]:.3e}, rate_rho_L2 ≈ {e1:.2f},rate_rho_H1 ≈ {e2:.2f}")
print(f"h={hs[j-1]:.3e} → h={hs[j]:.3e}, rate_u_L2 ≈ {r1:.2f},rate_u_H1 ≈ {r2:.2f}")
print(f"\nTotal runtime: {elapsed:.2f} seconds")
Why does the error convergence order of ρ only reach 2nd order, and not 3rd order? Looking forward to your reply！