Hello,I’m a fenicsx beginner and I’m trying to solve the maxwell equation in 2D.
I use C-N format for time discretization and finite element for space discretization.
I use n1curl elements with 2 degrees of freedom for the electric field E and CG elements with 1 degrees of freedom for the magnetic field H.
In theory, the convergence rate of both E and H should be 2, but the convergence rate of E is only 1.5, and the convergence rate of H is 2.
However, when I set H degrees of freedom to 2, E converges to 2, but H is only 2.
What’s wrong with my code? Is there any way that H and E can achieve the same convergence rate and conform to the theory?
Looking forward to your reply
Equation
Code
Here is my code
from mpi4py import MPI
import os
import math
import numpy as np
from petsc4py import PETSc
import ufl
import scipy.special as sps
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot,fem
from dolfinx.fem import Function, functionspace
from dolfinx.fem.petsc import LinearProblem,assemble_matrix,assemble_vector
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square,exterior_facet_indices
from ufl import dx, grad, inner,curl,MixedFunctionSpace
Nx_list = [24,32,40,48]
E_errlist,H_errlist = [],[]
alphajm = 0.5
gamma_2_alpha,gamma_3_alpha,gamma_4_alpha = sps.gamma(2 - alphajm),sps.gamma(3 - alphajm),sps.gamma(4 - alphajm)
T = 1
def E_expression(x,t):
values = np.zeros((2, x.shape[1]), dtype=np.float64)
values[0,:] = -np.sin(np.pi*x[1])*np.cos(np.pi*x[0])*(2*t**(2-alphajm)/gamma_3_alpha+t**2)
values[1,:] = np.cos(np.pi*x[1])*np.sin(np.pi*x[0])*(2*t**(2-alphajm)/gamma_3_alpha+t**2)
return values
def H_expression(x,t):
return -(2*t**(3-alphajm)/gamma_4_alpha+1/3*t**3)*2*np.pi*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])
def f_expression(x,t):
values = np.zeros((2, x.shape[1]), dtype=np.float64)
values[0,:] = (2*t**(1-alphajm)/gamma_2_alpha+2*t+2*np.pi**2*(2*t**(3-alphajm)/gamma_4_alpha+1/3*t**3))*(-np.cos(np.pi*x[0])*np.sin(np.pi*x[1]))
values[1,:] = (2*t**(1-alphajm)/gamma_2_alpha+2*t+2*np.pi**2*(2*t**(3-alphajm)/gamma_4_alpha+1/3*t**3))*(np.sin(np.pi*x[0])*np.cos(np.pi*x[1]))
return values
def get_b(t_temp,w_n):
E_n,H_n = w_n.split()
t1 = t_temp+dt/2
f_n = fem.Function(V3)
f_n.interpolate(lambda x: f_expression(x,t1))
b = (inner(E_n,v1)/dt+inner(curl(H_n),v1)/2+inner(f_n,v1))*dx+\
(inner(H_n,v2)/dt-inner(E_n,curl(v2))/2)*dx
return b
for Nx in Nx_list:
Ny=Nx
print(f"Nx : {Nx}")
dt=(1/Nx)
Mn = int(T/dt)
msh = create_unit_square(MPI.COMM_WORLD, Nx, Ny, CellType.quadrilateral)#CellType.quadrilateral
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
dx = ufl.Measure("dx", domain=msh, metadata={"quadrature_degree": 4})
curl1 = element("N1curl", msh.basix_cell(),2, dtype=default_real_type)
curlhigh = element("N1curl", msh.basix_cell(),4, dtype=default_real_type)
P1 = element("CG", msh.basix_cell(), 1, dtype=default_real_type)
Phigh = element("CG", msh.basix_cell(), 3, dtype=default_real_type)
mixV = functionspace(msh, mixed_element([curl1, P1]))
V1,V2,V3,V4 = functionspace(msh, curl1),functionspace(msh, P1),functionspace(msh, curlhigh),functionspace(msh, Phigh)
v1,v2 = ufl.TestFunctions(mixV)
E_u,H_u = ufl.TrialFunctions(mixV)
x = ufl.SpatialCoordinate(msh)
w_n = Function(mixV)
w_n.sub(0).interpolate(lambda x: E_expression(x, 0))
w_n.sub(1).interpolate(lambda x: H_expression(x, 0))
b = get_b(0,w_n)
a = (inner(E_u,v1)/dt-inner(curl(H_u),v1)/2)*dx+\
(inner(H_u,v2)/dt+inner(E_u,curl(v2))/2)*dx
for t in range(1,Mn+1):
t_temp = t*dt
#print(f't_temp: {t_temp }')
problem = LinearProblem(a, b, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
w_h = problem.solve()
E_0,H_0 = w_h.split()
b = get_b(t_temp,w_h)
solution = Function(mixV)
solution.sub(0).interpolate(lambda x: E_expression(x, t_temp))
solution.sub(1).interpolate(lambda x: H_expression(x, t_temp))
E_ex,H_ex = solution.split()
E_local = fem.assemble_scalar(fem.form(ufl.inner(E_0 - E_ex,E_0 - E_ex) *dx))
H_local = fem.assemble_scalar(fem.form(ufl.inner(H_0 - H_ex,H_0 - H_ex) *dx))
E_L2 = np.sqrt(msh.comm.allreduce(E_local, op=MPI.SUM))
H_L2 = np.sqrt(msh.comm.allreduce(H_local, op=MPI.SUM))
E_errlist.append(E_L2)
H_errlist.append(H_L2)
print(f"E_L2 : {E_L2:.4e}")
print(f"H_L2 : {H_L2:.4e}")
for i in range(len(Nx_list)-1):
orderE = np.log(E_errlist[i]/E_errlist[i+1])/np.log(Nx_list[i+1]/Nx_list[i])
orderH = np.log(H_errlist[i]/H_errlist[i+1])/np.log(Nx_list[i+1]/Nx_list[i])
print(f"orderE : {orderE:.4e}",f"orderH : {orderH:.4e}")