# Error analysis of decoupling the calculation of the MHD equation

**URL:** <https://fenicsproject.discourse.group/t/error-analysis-of-decoupling-the-calculation-of-the-mhd-equation/18086>\
**Category:** General\
**Tags:** mesh\
**Created:** [August 8, 2025, 8:21am UTC](https://fenicsproject.discourse.group/t/error-analysis-of-decoupling-the-calculation-of-the-mhd-equation/18086 "2025-08-08T08:21:29Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Dainy](https://avatars.discourse-cdn.com/v4/letter/d/b9bd4f/32.png) [@Dainy](https://fenicsproject.discourse.group/u/Dainy)\
**Post date:** [August 8, 2025, 8:21am UTC](https://fenicsproject.discourse.group/t/error-analysis-of-decoupling-the-calculation-of-the-mhd-equation/18086/1 "2025-08-08T08:21:30Z")

</div>

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

```auto
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:

```auto
(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?

```auto
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!

---

<div class="post-metadata">

**Author:** ![thingway](https://avatars.discourse-cdn.com/v4/letter/t/a5b964/32.png) [@thingway](https://fenicsproject.discourse.group/u/thingway)\
**Post date:** [November 11, 2025, 7:54am UTC](https://fenicsproject.discourse.group/t/error-analysis-of-decoupling-the-calculation-of-the-mhd-equation/18086/3 "2025-11-11T07:54:23Z")

</div>

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` .