Convergence rate of Maxwell simulation using custom L2 norm and cell-center sampling

mathematics
1

I am currently implementing Maxwell’s equations with a Cohen Monk perfectly matched layer in a two-dimensional time domain using a rectangular grid. The integration points for error calculation are only taken at the center of each element. Theoretically, the convergence order should be 2, but I consistently obtain an order of 1. I am not sure where the issue lies.Here is my code：

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import basix

from basix import CellType
from pathlib import Path
from dolfinx import fem, default_real_type, mesh
from dolfinx.fem import (
    Constant,
    Function,
    extract_function_spaces,
    functionspace,
    assemble_scalar,
    dirichletbc,
    form,
    Expression,
    locate_dofs_geometrical,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import (
    assemble_matrix,
    assemble_vector,
    apply_lifting,
    create_vector,
    set_bc,
)
from dolfinx.io import VTXWriter,XDMFFile
from dolfinx.mesh import create_rectangle, compute_midpoints, locate_entities_boundary, locate_entities, CellType, create_unit_square
from dolfinx.plot import vtk_mesh
from basix.ufl import element,  quadrature_element
from ufl import (
    FacetNormal,
    Identity,
    TestFunction,
    TrialFunction,
    div,
    curl,
    dot,
    ds,
    dx,
    inner,
    lhs,
    nabla_grad,
    rhs,
    sym,
    interpolate,
    SpatialCoordinate,
    as_vector,
    Measure,
    exp,
    sin,
    cos,
    pi,
    as_matrix
)
def l2_error(solution, exact_expr, space_type="E"):
    # Dx = dx(metadata={'quadrature_degree':1})
    if space_type == "E":
        # 计算误差表达式
        error_expr = form(inner(solution - exact_expr, solution - exact_expr) * dx)
    else:
        # 对于标量磁场
        error_expr = form((solution - exact_expr) * (solution - exact_expr) * dx)
    
    error = np.sqrt(MPI.COMM_WORLD.allreduce(assemble_scalar(error_expr), op=MPI.SUM))
    
    
    return error


def maxwell_simulation(N, dt, T=1.0, output_dir="results"):
    output_folder = Path(output_dir) / f"N{N}"
    output_folder.mkdir(exist_ok=True, parents=True)

    ax0 = 0;ay0 = 0
    bx0 = 1;by0 = 1

    Nx , Ny=N, N
    points = np.array([[ax0,ay0],[bx0,by0]], dtype=np.float64)
    n = [Nx, Ny]

    mesh = create_rectangle(MPI.COMM_WORLD, points, n, CellType.quadrilateral)
    V = element("N1curl", mesh.basix_cell(), 1)
    U = element("DG", mesh.basix_cell(), 0)
    Vs = functionspace(mesh, V)
    Us = functionspace(mesh, U)
    V_lagrange = functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
    U_lagrange = functionspace(mesh, ("Lagrange", 1))
    x = SpatialCoordinate(mesh)

    # 使用 ufl 定义表达式
    def curl_2d_E(E):
        return ufl.Dx(E[1], 0) - ufl.Dx(E[0], 1)

    # 精确解 - 使用 ufl 表达式
    def E_exact_expr(x, t):
        return as_vector([
        ufl.exp(-ufl.pi * t) * ufl.cos(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]),
        -ufl.exp(-ufl.pi * t) * ufl.sin(ufl.pi * x[0]) * ufl.cos(ufl.pi * x[1])
    ])
    # 用于插值的 numpy 函数（需要传入时间参数）
    def E_exact_np(x, t):
        return np.array([
            np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]),
            -np.exp(-np.pi * t) * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1])
        ])
    

    def H_exact_expr(x, t):
        return ufl.exp(-ufl.pi * t) * ufl.cos(ufl.pi * x[0]) * ufl.cos(ufl.pi * x[1])
    def H_exact_np(x, t):
        return np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.cos(np.pi * x[1])
    

    def E_star_exact_expr(x, t):
        return as_vector([ufl.exp(-pi * t) * cos(pi * x[0]) * sin(pi * x[1]),
                          -ufl.exp(-pi * t) * sin(pi * x[0]) * cos(pi * x[1])])
    def E_star_exact_np(x, t):
        return np.array([
            np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]),
            -np.exp(-np.pi * t) * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1])
        ])
    
    
    def K_exact_expr(x, t):
        return -(1 / ufl.pi) * ufl.exp(-ufl.pi * t) * ufl.cos(ufl.pi * x[0]) * ufl.cos(ufl.pi * x[1])
    def K_exact_np(x, t):
        return -(1 / np.pi) * np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.cos(np.pi * x[1])
    
    
    def H_star_exact_np(x_val, t_val):
        return np.exp(-np.pi * t_val) * np.cos(np.pi * x_val[0]) * np.cos(np.pi * x_val[1])
    

    def f_exact_np(x, t):
        return np.array([
            np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * 
            ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2),
            np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * 
            ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2)
        ])
    
    # 源项 hat_f
    def f_hat_expr(x,t):
        return as_vector([
            (ufl.exp(-pi * t) - 1) * ufl.cos(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) * 
            ((ufl.sin(ufl.pi * x[1]))**2 - (ufl.sin(pi * x[0]))**2),
            (ufl.exp(-pi * t) - 1) * sin(pi * x[0]) * cos(pi * x[1]) * 
            ((ufl.sin(pi * x[1]))**2 - (ufl.sin(ufl.pi * x[0]))**2)
        ])
    
    def f_hat_np(x, t):
        return np.array([
            (np.exp(-np.pi * t) - 1) * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * 
            ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2),
            (np.exp(-np.pi * t) - 1) * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * 
            ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2)
        ])
    
    # 源项 g_z
    g_exact_expr = cos(pi * x[0]) * cos(pi * x[1]) * ((sin(pi * x[0]))**2 + (sin(pi * x[1]))**2)
    
    def g_exact_np(x, t):
        return np.cos(np.pi * x[0]) * np.cos(np.pi * x[1]) * ((np.sin(np.pi * x[0]))**2 + (np.sin(np.pi * x[1]))**2)
    
    # 源项 hat_g_z
    def g_hat_expr(t):
        return (ufl.exp(-pi * t) * cos(pi * x[0]) * cos(pi * x[1]) * 
                (((sin(pi * x[0]))**2 + (sin(pi * x[1]))**2) - 
                 (1 / pi) * (sin(pi * x[0]))**2 * sin(pi * x[1])**2) -
                cos(pi * x[0]) * cos(pi * x[1]) * ((sin(pi * x[0]))**2 + (sin(pi * x[1]))**2))
    
    def g_hat_np(x, t):
        return (np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.cos(np.pi * x[1]) * 
                (((np.sin(np.pi * x[0]))**2 + (np.sin(np.pi * x[1]))**2) - 
                 (1 / np.pi) * (np.sin(np.pi * x[0]))**2 * (np.sin(np.pi * x[1]))**2) -
                np.cos(np.pi * x[0]) * np.cos(np.pi * x[1]) * 
                ((np.sin(np.pi * x[0]))**2 + (np.sin(np.pi * x[1]))**2))
    
    # 源项 g_z^*
    def g_star_expr(t):
        return -3 * pi * ufl.exp(-pi * t) * cos(pi * x[0]) * cos(pi * x[1])
    
    def g_star_np(x, t):
        return -3 * np.pi * np.exp(-np.pi * t) * np.cos(np.pi * x[0]) * np.cos(np.pi * x[1])
    
    # sigma 函数
    sigma_x_expr = (sin(pi * x[0]))**2
    sigma_y_expr = (sin(pi * x[1]))**2
    
    sigma_x = sigma_x_expr
    sigma_y = sigma_y_expr
    
    # 矩阵系数
    C_2d = as_matrix([[sigma_y, 0],
                      [0, sigma_x]])
    G_2d = as_matrix([[sigma_x, 0],
                      [0, sigma_y]])
    C_1d = sigma_x + sigma_y
    D_1d = sigma_x * sigma_y
    
    # PEC边界条件
    def pec_boundary(x):
        tol = 1e-10
        return (np.abs(x[0] - 0) < tol) | (np.abs(x[0] - 1) < tol) | \
            (np.abs(x[1] - 0) < tol) | (np.abs(x[1] - 1) < tol)

    boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, pec_boundary)
    boundary_dofs_E = locate_dofs_topological(Vs, mesh.topology.dim - 1, boundary_facets)
    Ebc = dirichletbc(Function(Vs), boundary_dofs_E)
    Estarbc = dirichletbc(Function(Vs), boundary_dofs_E)

    
    # 初始值
    eps0 = 1.0  # ε0
    mu0 = 1.0   # μ0
    
    E_n = Function(Vs)
    E_n.interpolate(lambda x: E_exact_np(x, 0))
    
    E_starn = Function(Vs)
    E_starn.interpolate(lambda x: E_star_exact_np(x, 0))
    
    H_nhalf = Function(Us)
    H_nhalf.interpolate(lambda x: H_exact_np(x, 0.5 * dt))
    
    H_starn = Function(Us)
    H_starn.interpolate(lambda x: H_star_exact_np(x, 0.5 * dt))
    
    K_nhalf = Function(Us)
    K_nhalf.interpolate(lambda x: K_exact_np(x, 0.5 * dt))
    
    # 源项函数
    f_n = Function(Vs)
    f_n.interpolate(lambda x: f_exact_np(x, 0))
    
    f_hat = Function(Vs)
    f_hat.interpolate(lambda x: f_hat_np(x, 0.5 * dt))
    
    g_n = Function(Us)
    g_n.interpolate(lambda x: g_exact_np(x, 0))
    
    g_star = Function(Us)
    g_star.interpolate(lambda x: g_star_np(x, dt))
    
    g_hat = Function(Us)
    g_hat.interpolate(lambda x: g_hat_np(x, dt))

    #变分形式
    E_star = TrialFunction(Vs)
    psi = TestFunction(Vs)
    a1 = form(ufl.inner(E_star, psi) * dx )
    L1 = form(ufl.inner(E_starn, psi) * dx 
            + dt * ufl.inner(H_nhalf, curl(psi)) * dx)
    A1 = assemble_matrix(a1)
    A1.assemble()
    b1 = create_vector(Vs)
    E_starsolution = Function(Vs)

    E = TrialFunction(Vs)
    psi_hat = TestFunction(Vs)
    a2 = form(ufl.inner(E, psi_hat) * dx + 0.5 * dt * ufl.inner(C_2d * E, psi_hat) * dx )
    L2 = form(ufl.inner(E_n, psi_hat) * dx 
            - 0.5 * dt * ufl.inner(C_2d * E_n, psi_hat) * dx 
            + ufl.inner(E_starsolution, psi_hat) * dx 
            - ufl.inner(E_starn, psi_hat) * dx 
            + 0.5 * dt * ufl.inner(G_2d * E_starsolution, psi_hat) * dx 
            + 0.5 * dt * ufl.inner(G_2d * E_starn, psi_hat) * dx
            + dt * ufl.inner(f_n, psi_hat) * dx
            + dt * ufl.inner(f_hat, psi_hat) * dx)
    A2 = assemble_matrix(a2)
    A2.assemble()
    b2 = create_vector(Vs)
    Esolution = Function(Vs)

    H_star = TrialFunction(Us)
    phi = TestFunction(Us)
    a3 = form(ufl.inner(H_star, phi) * dx )
    L3 = form(ufl.inner(H_starn, phi) * dx
            - dt * ufl.inner(curl(Esolution), phi) * dx
            + dt * ufl.inner(g_star, phi) * dx)
    A3 = assemble_matrix(a3)
    A3.assemble()
    b3 = create_vector(Us)
    H_starsolution = Function(Us)

    H = TrialFunction(Us)
    phi_hat = TestFunction(Us)
    a4 = form(ufl.inner(H, phi_hat) * dx 
            + 0.5 * dt * ufl.inner(C_1d * H, phi_hat) * dx
            + 0.25 * (dt**2) * ufl.inner(D_1d * H, phi_hat) * dx)
    L4 = form(ufl.inner(H_nhalf, phi_hat) * dx 
            - 0.5 * dt * ufl.inner(C_1d * H_nhalf, phi_hat) * dx
            - dt * ufl.inner(D_1d * K_nhalf, phi_hat) * dx
            - 0.25 * (dt**2) * ufl.inner(D_1d * H_nhalf, phi_hat) * dx
            + ufl.inner(H_starsolution,phi_hat) * dx
            - ufl.inner(H_starn, phi_hat) * dx
            + dt * ufl.inner(g_n, phi_hat) * dx
            + dt * ufl.inner(g_hat, phi_hat) * dx)
    A4 = assemble_matrix(a4)
    A4.assemble()
    b4 = create_vector(Us)
    Hsolution = Function(Us)

    K = TrialFunction(Us)
    phi_wide = TestFunction(Us)
    a5 = form(ufl.inner(D_1d * K, phi_wide) * dx)
    L5 = form(ufl.inner(D_1d * K_nhalf, phi_wide) * dx
        + 0.5 * dt * ufl.inner(D_1d * Hsolution, phi_wide) * dx
        + 0.5 * dt * ufl.inner(D_1d * H_nhalf, phi_wide) * dx)
    A5 = assemble_matrix(a5)
    A5.assemble()
    b5 = create_vector(Us)
    Ksolution = Function(Us)

    solver1 = PETSc.KSP().create(mesh.comm)
    solver1.setOperators(A1)
    solver1.setType(PETSc.KSP.Type.PREONLY)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.LU)

    solver2 = PETSc.KSP().create(mesh.comm)
    solver2.setOperators(A2)
    solver2.setType(PETSc.KSP.Type.PREONLY)
    pc2 = solver2.getPC()
    pc2.setType(PETSc.PC.Type.LU)

    solver3 = PETSc.KSP().create(mesh.comm)
    solver3.setOperators(A3)
    solver3.setType(PETSc.KSP.Type.PREONLY)
    pc3 = solver3.getPC()
    pc3.setType(PETSc.PC.Type.LU)

    solver4 = PETSc.KSP().create(mesh.comm)
    solver4.setOperators(A4)
    solver4.setType(PETSc.KSP.Type.PREONLY)
    pc4 = solver4.getPC()
    pc4.setType(PETSc.PC.Type.LU)

    solver5 = PETSc.KSP().create(mesh.comm)
    solver5.setOperators(A5)
    solver5.setType(PETSc.KSP.Type.PREONLY)
    pc5 = solver5.getPC()
    pc5.setType(PETSc.PC.Type.LU)
 
    
    E_for_output = Function(V_lagrange)
    H_for_output = Function(U_lagrange)
    folder = output_folder
    # 写入网格文件
    with XDMFFile(mesh.comm, folder / "mesh.xdmf", "w") as xdmf_file:
        xdmf_file.write_mesh(mesh)

    E_for_output = Function(V_lagrange)

    xdmf_E = XDMFFile(mesh.comm, folder / "Esolution_over_time.xdmf", "w")
    xdmf_H = XDMFFile(mesh.comm, folder / "Hsolution_over_time.xdmf", "w")
    xdmf_E.write_mesh(mesh)
    xdmf_H.write_mesh(mesh)
        
    t_E = 0
    t_half = 0.5 * dt
    Nt = 1000

    for n in range(Nt):
        t_E_new = t_E + dt
        t_next_half = t_half + dt

        with b1.localForm() as loc_1:
            loc_1.set(0)
            assemble_vector(b1, L1)
            apply_lifting(b1, [a1], [[Estarbc]])
            b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
            set_bc(b1, [Estarbc])
            solver1.solve(b1, E_starsolution.x.petsc_vec)
            E_starsolution.x.scatter_forward()

        f_hat.interpolate(lambda x: np.array([
            (np.exp(-pi * t_half) - 1) * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2),
            (np.exp(-pi * t_half) - 1) * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * ((np.sin(np.pi * x[1]))**2 - (np.sin(np.pi * x[0]))**2)
            ]))

        with b2.localForm() as loc_2:
            loc_2.set(0)
            assemble_vector(b2, L2)
            apply_lifting(b2, [a2], [[Ebc]])
            b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
            set_bc(b2, [Ebc])
            solver2.solve(b2, Esolution.x.petsc_vec)
            Esolution.x.scatter_forward()

        g_star.interpolate(lambda x: g_star_np(x, t_E_new))

        with b3.localForm() as loc_3:
            loc_3.set(0)
        assemble_vector(b3, L3)
        apply_lifting(b3, [a3], [[]])
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b3, [])
        solver3.solve(b3, H_starsolution.x.petsc_vec)
        H_starsolution.x.scatter_forward()

        g_hat.interpolate(lambda x: g_hat_np(x, t_E_new))

        with b4.localForm() as loc_4:
            loc_4.set(0)
        assemble_vector(b4, L4)
        apply_lifting(b4, [a4], [[]])
        b4.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b4, [])
        solver4.solve(b4, Hsolution.x.petsc_vec)
        Hsolution.x.scatter_forward()

        with b5.localForm() as loc_5:
            loc_5.set(0)
        assemble_vector(b5, L5)
        b5.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b5, [])
        solver5.solve(b5, Ksolution.x.petsc_vec)
        Ksolution.x.scatter_forward()

        # 每100步或最后一步保存XDMF文件
        if n % 100 == 0 or n == Nt - 1:
            E_for_output.interpolate(Esolution)
            H_for_output.interpolate(Hsolution)
            E_for_output.name = "E_field"
            # Hsolution.name = "H_field"
            H_for_output.name = "H_field"

            # Hsolution.name = "H_field"
            xdmf_E.write_function(E_for_output, t_E_new)
            xdmf_H.write_function(H_for_output, t_next_half)
            # vtk_E.write_function(E_for_output, t_E_new)
            # vtk_H.write_function(Hsolution, t_next_half)
            print(f"  -> Data written at t={t_E_new:.6e} s (={t_E_new*1e9:.3f} ns)")
       
        H_nhalf.x.array[:] = Hsolution.x.array[:]
        E_n.x.array[:] = Esolution.x.array[:]
        H_starn.x.array[:] = H_starsolution.x.array[:]
        E_starn.x.array[:] = E_starsolution.x.array[:]
        K_nhalf.x.array[:] = Ksolution.x.array[:]
        E_n.x.scatter_forward()
        H_nhalf.x.scatter_forward()
        H_starn.x.scatter_forward()
        E_starn.x.scatter_forward()
        K_nhalf.x.scatter_forward()

        t_E = t_E_new
        t_half = t_next_half

    

    # 计算最终时刻误差
    E_expr_final = E_exact_expr(x, t_E)
    H_expr_final = H_exact_expr(x, t_half)
    
    final_E_error = l2_error(Esolution, E_expr_final, space_type="E")
    final_H_error = l2_error(Hsolution, H_expr_final, space_type="H")

    xdmf_E.close()
    xdmf_H.close()
    b1.destroy()
    b2.destroy()
    b3.destroy()
    b4.destroy()
    b5.destroy()
    solver1.destroy()
    solver2.destroy()
    solver3.destroy()
    solver4.destroy()
    
    

    h = 1.0/N  # 网格尺寸
    return {
        'N': N,
        'h': h,
        'dt': dt,
        'final_time_E': t_E,
        'final_time_H': t_half,
        'error_E_abs': final_E_error,
        'error_H_abs': final_H_error,
        'num_time_steps': Nt
    }

def calculate_convergence_rates(results):
    
    print("空间收敛性分析")
   
    # 提取数据
    N_values = [r['N'] for r in results]
    h_values = [r['h'] for r in results]
    # dt_values = [r['dt'] for r in results]
    E_errors = [r['error_E_abs'] for r in results]
    H_errors = [r['error_H_abs'] for r in results]
    
    print(f"{'N':>6} {'h':>10} {'E_err(绝对)':>15} {'H_err(绝对)':>15} ")
    print("-" * 85)
    
    for i in range(len(results)):
        print(f"{N_values[i]:6d} {h_values[i]:10.6f} "
              f"{E_errors[i]:15.6e} {H_errors[i]:15.6e} ")#{dt_values[i]:12.6e}
    
    # 计算收敛阶
    print("\n收敛阶计算:")
    print("-" * 50)
    
    convergence_rates = []
    for i in range(1, len(h_values)):
        rate_E = np.log(E_errors[i-1] / E_errors[i]) / np.log(h_values[i-1] / h_values[i])
        rate_H = np.log(H_errors[i-1] / H_errors[i]) / np.log(h_values[i-1] / h_values[i])
        print(f"h: {h_values[i-1]:.4f} -> {h_values[i]:.4f} (细化比: {h_values[i-1]/h_values[i]:.2f})")
        print(f"  E场绝对误差收敛阶: {rate_E:.4f}")
        print(f"  H场绝对误差收敛阶: {rate_H:.4f}")
    
    return {
        'N_values': N_values,
        'h_values': h_values,
        'E_errors': E_errors,
        'H_errors': H_errors,
        'convergence_rates': convergence_rates
    }

def spatial_convergence_test():
 
    print("="*70)
    print("开始空间收敛性测试")
    print("="*70)
    
    
    T = 0.01   # 总时间
    dt = 1e-5# 固定时间步长
    N_values = [10, 20, 40, 80, 160]
    results = []
    
    for N in N_values:
        h = 1/N
        # dt = 0.01 * h
        print(f"\n运行 N={N}, h={1.0/N:.4f}, dt={dt:.6e}, T={T}")
        print("-" * 40)
        
        result = maxwell_simulation(
            N=N,
            dt=dt,
            T=T,
            output_dir=f"cr2"
        )
        
        results.append(result)
        
        print(f"完成! 最终误差:")
        print(f"  E场绝对误差: {result['error_E_abs']:.6e}")
        print(f"  H场绝对误差: {result['error_H_abs']:.6e}")
    
    # 计算收敛阶
    convergence_data = calculate_convergence_rates(results)
    return {
        'results': results,
        'convergence_data': convergence_data
    }

if __name__ == "__main__":
    results = spatial_convergence_test()

Here are my results.


Estimated Convergence Rates
Refinement	E field order	H field order
0.1000 → 0.0500	1.4940	0.5751
0.0500 → 0.0250	1.3871	0.6573
0.0250 → 0.0125	1.0953	0.9593
0.0125 → 0.0063	0.9845	1.2414
2

This code is simply to lengthy for me to parse in detail.
Would you mind:

  1. Reducing it to solving just one of the four equation above, and use manufactured terms for all other terms (that in your code stems from solving the remaining 3 PDEs).
  2. Ensure that your finite element pair is stable for the discretization.
  3. It is unclear where your custom L2-norm enters, as there is no such thing in the code (that I can see at a first glance). If you mean by setting H in DG-0, that is not changing the norm (per say). If you want custom interpolation in the norm, you should pass the custom integration point (or set a low order quadrature degree) in the integration measure dx.
3

Sorry, I didn’t understand your meaning. The two-dimensional time-domain Maxwell equations with Cohen Monk PML consist of five equations. At the end of the time loop, I called the defined L² norm to calculate the errors of the electric and magnetic fields at the final time for convergence order analysis. This is my minimal code：

def maxwell_simulation(N, dt, T=0.01):
    domain = create_rectangle(MPI.COMM_WORLD, [[0,0], [1,1]], [N, N], CellType.quadrilateral)
    x = SpatialCoordinate(domain)
    
    V = functionspace(domain, ("N1curl", 1))
    U = functionspace(domain, ("DG", 0))
    
    sigma_x = (sin(pi * x[0]))**2
    sigma_y = (sin(pi * x[1]))**2
    C_2d = as_matrix([[sigma_y, 0], [0, sigma_x]])
    G_2d = as_matrix([[sigma_x, 0], [0, sigma_y]])
    C_1d = sigma_x + sigma_y
    D_1d = sigma_x * sigma_y
    
    def pec_boundary(x):
        tol = 1e-10
        return (np.abs(x[0]-0)<tol) | (np.abs(x[0]-1)<tol) | (np.abs(x[1]-0)<tol) | (np.abs(x[1]-1)<tol)
    
    boundary_facets = locate_entities_boundary(domain, domain.topology.dim-1, pec_boundary)
    boundary_dofs = fem.locate_dofs_topological(V, domain.topology.dim-1, boundary_facets)
    bc = dirichletbc(Function(V), boundary_dofs)
    
    E_star, E = TrialFunction(V), TrialFunction(V)
    psi, psi_hat = TestFunction(V), TestFunction(V)
    
    H_star, H = TrialFunction(U), TrialFunction(U)
    phi, phi_hat = TestFunction(U), TestFunction(U)
    
    K = TrialFunction(U)
    phi_wide = TestFunction(U)
    
    a1 = form(inner(E_star, psi) * dx)
    L1 = form(inner(Function(V), psi) * dx + dt * inner(Function(U), curl(psi)) * dx)
    
    a2 = form(inner(E, psi_hat) * dx + 0.5 * dt * inner(C_2d * E, psi_hat) * dx)
    L2 = form(inner(Function(V), psi_hat) * dx - 0.5 * dt * inner(C_2d * Function(V), psi_hat) * dx
            + inner(Function(V), psi_hat) * dx - inner(Function(V), psi_hat) * dx
            + 0.5 * dt * inner(G_2d * Function(V), psi_hat) * dx + 0.5 * dt * inner(G_2d * Function(V), psi_hat) * dx
            + dt * inner(Function(V), psi_hat) * dx + dt * inner(Function(V), psi_hat) * dx)
    
    a3 = form(inner(H_star, phi) * dx)
    L3 = form(inner(Function(U), phi) * dx - dt * inner(curl(Function(V)), phi) * dx
            + dt * inner(Function(U), phi) * dx)
    
    a4 = form(inner(H, phi_hat) * dx + 0.5 * dt * inner(C_1d * H, phi_hat) * dx
            + 0.25 * (dt**2) * inner(D_1d * H, phi_hat) * dx)
    L4 = form(inner(Function(U), phi_hat) * dx - 0.5 * dt * inner(C_1d * Function(U), phi_hat) * dx
            - dt * inner(D_1d * Function(U), phi_hat) * dx - 0.25 * (dt**2) * inner(D_1d * Function(U), phi_hat) * dx
            + inner(Function(U), phi_hat) * dx - inner(Function(U), phi_hat) * dx
            + dt * inner(Function(U), phi_hat) * dx + dt * inner(Function(U), phi_hat) * dx)
    
    a5 = form(inner(D_1d * K, phi_wide) * dx)
    L5 = form(inner(D_1d * Function(U), phi_wide) * dx + 0.5 * dt * inner(D_1d * Function(U), phi_wide) * dx
            + 0.5 * dt * inner(D_1d * Function(U), phi_wide) * dx)
    
    A1, A2, A3, A4, A5 = [assemble_matrix(a) for a in [a1, a2, a3, a4, a5]]
    [A.assemble() for A in [A1, A2, A3, A4, A5]]
    
    b1, b2, b3, b4, b5 = [create_vector(V if i<2 else U) for i in range(5)]
    solutions = [Function(V if i<2 else U) for i in range(5)]
    
    solvers = []
    for A in [A1, A2, A3, A4, A5]:
        solver = PETSc.KSP().create(domain.comm)
        solver.setOperators(A)
        solver.setType(PETSc.KSP.Type.PREONLY)
        solver.getPC().setType(PETSc.PC.Type.LU)
        solvers.append(solver)
    
    Nt = int(T/dt)
    for n in range(Nt):
        for i, (A, b, L, solver, sol, bc_list) in enumerate([
            (A1, b1, L1, solvers[0], solutions[0], [bc]),
            (A2, b2, L2, solvers[1], solutions[1], [bc]),
            (A3, b3, L3, solvers[2], solutions[2], []),
            (A4, b4, L4, solvers[3], solutions[3], []),
            (A5, b5, L5, solvers[4], solutions[4], [])
        ]):
            with b.localForm() as loc:
                loc.set(0)
            assemble_vector(b, L)
            apply_lifting(b, [A], [bc_list])
            b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
            set_bc(b, bc_list)
            solver.solve(b, sol.x.petsc_vec)
            sol.x.scatter_forward()
    
    return {'N': N, 'h': 1.0/N, 'dt': dt}

def spatial_convergence_test():
    dt = 1e-5
    N_values = [10, 20, 40, 80, 160]
    results = []
    
    for N in N_values:
        result = maxwell_simulation(N=N, dt=dt)
        results.append(result)
    
    return results

if __name__ == "__main__":
    spatial_convergence_test()

And I’ve tried custom integration rules,

dx_center = ufl.Measure(
#         "dx",
#         domain=mesh,
#         metadata={
#             "quadrature_scheme": "custom",
#             "quadrature_points": points,
#             "quadrature_weights": weights,
#         }
#     )

but the results appear to be the same as when I simply specify an integration precision of 1 like below:

Dx = dx(metadata={'quadrature_degree':1})

The equations I’m solving are as follows：

\left(\sigma_{1}\sigma_{2}\,\frac{K_{h}^{n+\frac{3}{2}}-K_{h}^{n+\frac{1}{2}}}{\tau},\widetilde{\Phi}_{h}\right) =\left(\sigma_{1}\sigma_{2}\,\frac{H_{h}^{n+\frac{3}{2}}+H_{h}^{n+\frac{1}{2}}}{2},\widetilde{\Phi}_{h}\right), \quad\forall\widetilde{\Phi}_{h}\in U_{h}, \tag{1}
\left(\epsilon_{0}\frac{E_{h}^{*n+1}-E_{h}^{*n}}{\tau},\Psi_{h}\right) =\left(H_{h}^{n+\frac{1}{2}},\nabla\times\Psi_{h}\right), \quad\forall\Psi_{h}\in V_{h}^{0}, \tag{2}
\left(\mu_{0}\frac{H_{h}^{*n+\frac{3}{2}}-H_{h}^{*n+\frac{1}{2}}}{\tau},\Phi_{h}\right) =-\left(\nabla\times E_{h}^{n+1},\Phi_{h}\right) +\left(g_{z}^{*}\left(t_{n+1}\right),\Phi_{h}\right), \quad\forall\Phi_{h}\in U_{h}, \tag{3}
\left(\frac{E_{h}^{n+1}-E_{h}^{n}}{\tau},\widehat{\Psi}_{h}\right) +\left(C_{2d}\,\frac{E_{h}^{n+1}+E_{h}^{n}}{2},\widehat{\Psi}_{h}\right)\\ = \left(\frac{E_{h}^{*n+1}-E_{h}^{*n}}{\tau},\widehat{\Psi}_{h}\right) +\left(G_{2d}\frac{E_{h}^{*n+1}+E_{h}^{*n}}{2},\widehat{\Psi}_{h}\right) +\left((f+\hat{f})(t_{n+\frac{1}{2}}),\widehat{\Psi}_{h}\right), \quad\forall\widehat{\Psi}_{h}\in V_{h}^{0}, \tag{4}
\left(\frac{H_{h}^{n+\frac{3}{2}}-H_{h}^{n+\frac{1}{2}}}{\tau},\widehat{\Phi}_{h}\right) +\left((\sigma_{1}+\sigma_{2})\frac{H_{h}^{n+\frac{3}{2}}+H_{h}^{n+\frac{1}{2}}}{2},\widehat{\Phi}_{h}\right)+\left(\sigma_{1}\sigma_{2}\frac{K_{h}^{n+\frac{3}{2}}+K_{h}^{n+\frac{1}{2}}}{2},\widehat{\Phi}_{h}\right)\\ =\left(\frac{H_{h}^{*n+\frac{3}{2}}-H_{h}^{*n+\frac{1}{2}}}{\tau},\widehat{\Phi}_{h}\right) +\left((g_{z}+\hat{g}_{z})\left(t_{n+1}\right),\widehat{\Phi}_{h}\right), \quad\forall\widehat{\Phi}_{h}\in U_{h}. \tag{5}
4

My point is that since you have a manufactured solution, you can reduce say the first system

F(K^{n+\frac{3}{2}}, K^{n+\frac{1}{2}}, H^{n+\frac{3}{2}}, H^{n+\frac{1}{2}}, \tilde \phi_h)=0
to
F(K^{n+\frac{3}{2}}, K_{ex}^{n+\frac{1}{2}}, H_{ex}^{n+\frac{3}{2}}, H_{ex}^{n+\frac{1}{2}}, \tilde \phi_h)=0

where the _{ex} subscript indicates that you have inserted the manufactured/exact solution for K and H. If you solve that for a single time step, do you observe the expected convergence in K?

Once it converges for a single time step, check convergence for multiple steps.

Then, step by step you can replace the components of your segregated solver (splitting scheme) with each of these steps.

As I am not familiar with the particular splitting scheme that you are using to solve these equations, its hard for me to comment on the overall observed behavior.

5

Thank you for your reply. I will make another attempt to revise it.