Error analysis of decoupling the calculation of the MHD equation

General
1

Hi, Everyone！
Currently, I am decoupling and computing a 3D model：

E+ u\times B-(1/Rm)\curl B=g  in \Omega         E\times n=0    on \partial Omaga

Here, E represents the electric field, B represents the magnetic field, and u represents the velocity field. The discrete format of the equation is:

(E_h^n,F)+(u_h^n\times B_h^{n-1},F)-(1/Rm)*(B_h^n,\curl F)=(g,F)

I use the RT0 element to approximate B_h^n, the ND0(lowest-order edge element) to approximate E_h^n, and the vector P1 element to approximate u. And the B and u involved in the code are all calculated using the exact solution, but the convergence order I calculated was always incorrect. Can you help me see what’s wrong with the code?

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
import time 
from ufl import dx, grad, dot, lhs, rhs,inner,cross,curl,ds
from basix.ufl import mixed_element,element
from numpy import sin, cos, pi, exp
# 精确解
Rm=1
left =0.0
right = 1.0  
bottom =  0.0
top = 1.0     
near = 0.0
far = 1.0
def  ExactU(x,t):
    factor = exp(-t) + 1.0
    temp1=factor * sin(pi * x[0]) *sin(pi * x[1]) *sin(pi * x[2])
    temp2=np.zeros_like(x[0])
    temp3=factor * (x[2]**2 - x[2]) *sin(pi * x[0]) * sin(pi * x[1])
    return [temp1,temp2,temp3]


def ExactE(x,t):
    sint = sin(pi*t)+1.5
    temp1=sint*(x[1]**2-x[1])*sin(pi*x[2])
    temp2=sint*(x[0]**2-x[0])*sin(pi*x[2])
    temp3=sint*(x[0]**2-x[0])*sin(pi*x[1])
    return [temp1,temp2,temp3]




def ExactB(x,t):
    cost=(1/pi)*cos(pi*t)-1.5*t
    temp1=cost*pi*(x[0]**2-x[0])*(cos(pi*x[1])-cos(pi*x[2]))
    temp2=cost*(pi*(x[1]**2-x[1])*cos(pi*x[2])-(2*x[0]-1)*sin(pi*x[1]))
    temp3=cost*(2*x[0]-2*x[1])*sin(pi*x[2])
    return [temp1,temp2,temp3]



def rhs_expr(x, t):
      factor = np.exp(-t) + 1.0
      sint=sin(pi*t)+1.5
      cost=1/pi*cos(pi*t)-1.5*t

      E1=(x[1]**2-x[1])*sin(pi*x[2])
      E2=(x[0]**2-x[0])*sin(pi*x[2])
      E3=(x[0]**2-x[0])*sin(pi*x[1])

      u1= sin(pi * x[0]) * sin(pi * x[1]) *sin(pi * x[2])
      u2= np.zeros_like(x[0])
      u3= (x[2]**2 - x[2]) * sin(pi * x[0]) * sin(pi * x[1])

      B1=pi*(x[0]**2-x[0])*(cos(pi*x[1])-cos(pi*x[2]))
      B2=(pi*(x[1]**2-x[1])*cos(pi*x[2])-(2*x[0]-1)*sin(pi*x[1]))
      B3=(2*x[0]-2*x[1])*sin(pi*x[2])

      curlB1=-2*sin(pi*x[2])+pi**2*(x[1]**2-x[1])*sin(pi*x[2])
      curlB2=-2*sin(pi*x[2])+pi**2*(x[0]**2-x[0])*sin(pi*x[2])
      curlB3=-2*sin(pi*x[1])+pi**2*(x[0]**2-x[0])*sin(pi*x[1])

      temp1=sint*E1\
            +factor*cost*(u2*B3-u3*B2)\
            -(1/Rm)*cost*curlB1
            
      temp2=sint*E2\
            +factor*cost*(u3*B1-u1*B3)\
            -(1/Rm)*cost*curlB2
           
         #   
      temp3=sint*E3\
            +factor*cost*(u1*B2-u2*B1)\
           -(1/Rm)*cost*curlB3
      return [temp1, temp2,temp3]

       
# 固定终止时间
T = 1 / 2
# 结果记录
hs = []
errors_L2 = []
errors_max = []  
errors_H1 = []
errors_curl =[]
start_time = time.time()
for i in range(3):
    n = 2 ** (i + 2)               # 空间划分数 n = 2^{i+2}
    h = 1.0 / n
    num_steps = 2 ** (2 * i + 3)   # 时间划分数 tau = 2^{2i+4}
    dt = T / num_steps
    t = 0.0
# 创建 3D 网格（单位立方体）
    domain = mesh.create_unit_cube(MPI.COMM_WORLD, n, n, n, mesh.CellType.tetrahedron)
    x = ufl.SpatialCoordinate(domain)
    V = fem.functionspace(domain, ("Lagrange", 1))  # 线性 Lagrange 元
    Vv=fem.functionspace(domain, ("Lagrange", 1, (3, )))
#  磁场空间与电场空间的定义
    VB = fem.functionspace(domain, ("RT", 1))
    VE = fem.functionspace(domain, ("N1curl", 1)) 

    E_n=fem.Function(VE)
    u_n = fem.Function(Vv)

    B_old = fem.Function(VB)
    B_old.interpolate(lambda x: ExactB(x, t))

    B_n = fem.Function(VB)

    E_D = fem.Function(VE)
    E_D.interpolate(lambda x: ExactE(x, t)) 

    f_fun = fem.Function(VE)
    
# 施加 Dirichlet 边界条件（3D 边界）
    tdim = domain.topology.dim  
    fdim = tdim - 1  
    domain.topology.create_connectivity(fdim, tdim)
    boundary_facets = mesh.exterior_facet_indices(domain.topology)
    bc = fem.dirichletbc(E_D, fem.locate_dofs_topological(VE, fdim, boundary_facets))

    ns = ufl.FacetNormal(domain)
# 变分形式（3D 形式相同，但 dx 现在是 3D 积分）
    E, F = ufl.TrialFunction(VE), ufl.TestFunction(VE)
    G = inner(E,F)*dx + inner(cross(u_n,B_old),F)*dx- inner(f_fun,F)*dx \
        - (1/Rm)*inner(B_n,curl(F))*dx 
    
    a = fem.form(lhs(G))
    L = fem.form(rhs(G))

# 组装矩阵和向量
    A = assemble_matrix(a, bcs=[bc])
    A.assemble()
    b = create_vector(L)
    E_h =fem.Function(VE) 

# 定义求解器（直接法：LU 分解）
    solver = PETSc.KSP().create(domain.comm)
    solver.setOperators(A)
  #  solver.setType(PETSc.KSP.Type.PREONLY)
  #  solver.getPC().setType(PETSc.PC.Type.LU)

    solver.setType(PETSc.KSP.Type.BCGS)
    pc1 = solver.getPC()
    pc1.setType(PETSc.PC.Type.HYPRE)
    pc1.setHYPREType("boomeramg")


# 时间步进
    for n_step in range(num_steps):
      t += dt
      u_n.interpolate(lambda x:ExactU(x,t))
      B_n.interpolate(lambda x: ExactB(x, t))
      E_D.interpolate(lambda x: ExactE(x, t)) 
      bc = fem.dirichletbc(E_D, fem.locate_dofs_topological(VE, fdim, boundary_facets))
      f_fun.interpolate(lambda x: rhs_expr(x, t))
      # 更新矩阵
      A.zeroEntries()
      assemble_matrix(A, a, bcs=[bc])
    #  assemble_matrix(A, a, bcs=[bc])
      A.assemble()

      # 重新组装右端向量
      with b.localForm() as loc_b:
        loc_b.set(0)
      assemble_vector(b, L)

      # 施加边界条件
      apply_lifting(b, [a], [[bc]])
      b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
      set_bc(b, [])

      # 求解线性系统
      solver.solve(b, E_h.x.petsc_vec)
      E_h.x.scatter_forward()
      B_old.interpolate(lambda x: ExactB(x, t))
      
    
     

# 计算误差（在更高阶空间）
    V_dg = fem.functionspace(domain, ("Lagrange", 1, (3,)))
  #  V_dg = fem.functionspace(domain, ("N1curl", 2))
    E_ex = fem.Function(V_dg)
    E_ex.interpolate(lambda x: ExactE(x,t))

    E_ex1 = fem.Function(VE)
    E_ex1.interpolate(lambda x:ExactE(x,t))

# L2 误差
    error_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((E_h - E_ex)**2 * ufl.dx)), op=MPI.SUM))
    error_max= np.max(np.abs(E_ex1.x.array - E_h.x.array))
    error_curl = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form(curl(E_h - E_ex)**2*dx)), op=MPI.SUM))
    if MPI.COMM_WORLD.rank == 0:
        hs.append(h)
        errors_L2.append(error_L2)
        errors_max.append(error_max)
        errors_curl.append(error_curl)
        print(f"[i={i}] n={n}, num_steps={num_steps}, dt={dt:.3e}, h={h:.3e}")
        print(f"        → L2 error = {error_L2:.3e}, max error={error_max:.3e},curl error={error_curl:.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)):
        r1 = np.log(errors_L2[j-1]/errors_L2[j]) / np.log(hs[j-1]/hs[j])
        r2 = np.log(errors_max[j-1]/errors_max[j]) / np.log(hs[j-1]/hs[j])
        r3 = np.log(errors_curl[j-1]/errors_curl[j]) / np.log(hs[j-1]/hs[j])
        print(f"h={hs[j-1]:.3e} → h={hs[j]:.3e}, rate ≈ {r1:.2f},rate ≈ {r2:.2f},rate ≈ {r3:.2f}")

    print(f"\nTotal runtime: {elapsed:.2f} seconds")

Looking forward to your reply!

3

There are some advices

  1. In your code, V_dg is used to compute error_L2 and error_curl , but there are two different definitions for it. Based on my experience, it is reasonable to use a function space that is consistent with the one used for the numerical solution.
  2. In DefElement , if you intend to use the first-kind Nédélec element, the option in basic.ufl is N1E . You might consider replacing N1curl with N1E , although the precise differences between them are unclear to me.
  3. It appears that your code only solves for the electric field E numerically, while the velocity u_n and magnetic field B_n are obtained by interpolating the exact solution at certain times. I understand this is a sample code, and the full system includes equations for velocity and magnetic field. However, if the goal is to test the solver specifically for the electric field, it is acceptable to use the exact solutions for u and B by integrating them into the source term g .
  4. Your code uses a first-order Nédélec element. Theoretically, this predicts a first-order spatial convergence rate, which matches the expected temporal convergence rate (CR). Therefore, to observe a first-order convergence in the H(curl) norm, you should set the time step as Δt = h . However, your current code does not use it even Δt = h² because the final time is set to T = 1/2 .