Hi,
I’m using Fenics 2017.2.0 I want to test the convergence rate for the standard Maxwell’s equation in 2D case, where H is a vector and E is a scalar. The equations are given by:
\epsilon_0 \partial_t E_z = \nabla\times \mathbf{H} + f_1\\
\mu_0 \partial_t \mathbf{H} = -\nabla\times E_z +\mathbf{f_2}
The Crank-Nicolson scheme is given by following
\epsilon_0 (\frac{E^{n+1}_h-E^n_h}{\tau}, \phi_h) = (\nabla\times \frac{\mathbf{H}^{n+1}_h+\mathbf{H}^n_h}{2}, \phi_h) + (f_1^{n+\frac{1}{2}}, \phi_h) \\
\mu_0 (\frac{\mathbf{H}^{n+1}_h-\mathbf{H}^n_h}{\tau}, \mathbf{\psi_h}) =-(\frac{E^{n+1}_h+E^n_h}{2},\nabla\times \mathbf{\psi_h})+ (\mathbf{f_2}^{n+\frac{1}{2}}, \psi_h)
I got the right convergence rate for O(h^r), where r is the N2curl basis function order. But I can only observe O(\tau) instead of O(\tau^2) by the theoretical result. I can’t figure out why.
Really hope to get help from you. I attached my full code as follows:
# CRCN_MaxwellStandard
# Test for CN scheme with standard Maxwell's Equation
# Standard Maxwell's equation
# Want to observe O(dt^2 + h^r)
# I can see O(h^r) by fixing dt=h^r then let h vary
# I can't see O(dt^2) by two methods I wrote in the code but I can only see O(dt)
from copy import deepcopy
from fenics import *
from dolfin import *
from cmath import pi
import sympy as sym
from sympy import sin, cos, exp
import numpy as np
set_log_active(False)
# Use symbolic differentiation and integration to calculate J, K and f1, f2
x,y,t = sym.symbols('x[0], x[1], t')
omega =2
omega_t = 2
mu0=1.0; eps0=1.0;Cv=1.0
#Define the exact solution by symbolic tool
Hx = cos(omega*pi*x)*sin(omega*pi*y)*sin(omega_t*pi*t)*pow(t,1)
Hy = -sin(omega*pi*x)*cos(omega*pi*y)*sin(omega_t*pi*t)*pow(t,1)
E = cos(omega*pi*x)*cos(omega*pi*y)*sin(omega_t*pi*t)*pow(t,1)
# Define source term
f1 = eps0*sym.diff(E, t) - sym.diff(Hy, x) + sym.diff(Hx, y)
f2x = mu0* sym.diff(Hx, t) + sym.diff(E, y)
f2y = mu0* sym.diff(Hy, t) - sym.diff(E, x)
f1_code = sym.printing.ccode(f1)
f2x_code = sym.printing.ccode(f2x)
f2y_code = sym.printing.ccode(f2y)
E_code = sym.printing.ccode(E)
Hx_code = sym.printing.ccode(Hx)
Hy_code = sym.printing.ccode(Hy)
degree=4 # degree for exact solution expression
polyn =4 # degree for trail and test function coef expression
r=2 # order of edge element basis functions
Hsol = Expression((Hx_code,Hy_code), t=0, degree=degree)
Esol = Expression(E_code,t=0,degree=degree)
#vacuum domain
ax = ay =0.0
bx = by =1.0
harray = [] # mesh sizes
Error = [] # errors
Err = [] # errors
for i in range(1,6):
#test O(h^r)
# N = int(pow(2,i-1)*4);
# h=bx/N;harray.append(h)
# T = 1; dt = pow(h,r)
# Nt= int(T/dt)
# test O(dt^2), method1
N = int(pow(2,i-1)*4);
h=bx/N;harray.append(h)
T = 1; dt = h/4
Nt= int(T/dt)
print(Nt)
# test O(dt^2), method2
# N=64; h = bx/N
# harray.append(h)
# T = 1
# Nt = int(pow(2,i-1)*4)
# dt = T/Nt
# print(Nt)
mesh = RectangleMesh(Point(ax, ay), Point(bx, by), N, N)
# meshfile = File("CNMaxwellresult/mesh.pvd") #put the mesh data to a new folder
# meshfile << mesh
V = FiniteElement("N2curl", mesh.ufl_cell(), r)
U = FiniteElement("DG", mesh.ufl_cell(), r-1)
All_element = MixedElement([V, U])
All_space = FunctionSpace(mesh, All_element)
# define boundary condition
def boundary(x, on_boundary):
return on_boundary
Hbc = DirichletBC(All_space.sub(0), Constant((0.0,0.0)), boundary)
Ebc = DirichletBC(All_space.sub(1), Constant(0.0), boundary)
bcs = [Hbc,Ebc]
# Define initial condition, which is current step
solution00 = Constant((0.0, 0.0, 0.0))
solution0 = interpolate(solution00, All_space)
H, E = solution0.split(deepcopy=True)
# Define TrialFunction, which is next step
solution = TrialFunction(All_space)
Hn,En = split(solution)
# Define test functions
psi, phi= TestFunctions(All_space)
f1 = Expression(f1_code, t=0.5*dt ,degree=degree)
f2 = Expression((f2x_code,f2y_code),t=0.5*dt, degree=degree)
F = eps0*inner(En-E,phi)*dx - dt/2*inner(curl(Hn+H), phi)*dx - dt*inner(f1, phi)*dx\
+mu0*inner(Hn-H, psi)*dx + dt/2*inner(En+E,curl(psi))*dx - dt*inner(f2, psi)*dx
a,L = lhs(F), rhs(F)
# Define Function, which is unknown in the linear system need to be solved
solution_n = Function(All_space)
# Time-stepping
# vtkfile = File('CRM/Esolution.pvd')
# vtkfile1 = File('CRM/Eextsolution.pvd')
# vtkfile2 = File('CRM/Hsolution.pvd')
# vtkfile3 = File('CRM/Hextsolution.pvd')
t = 0
for n in range(Nt+1):
# print( 'time step is: %d (total steps: %d)' % (n, Nt))
# print("time", t)
solve(a == L, solution_n, bcs)
Hsolution, Esolution = solution_n.split(deepcopy=True)
solution0.assign(solution_n)
H,E = solution0.split(deepcopy=True)
t += dt
f1.t=t+0.5*dt
f2.t=t+0.5*dt
# vtkfile << Esolution
# vtkfile2 << Hsolution
print( 'mesh N: %d with r= %d and final time dt= %.6f' %(N,r,dt) )
print("final time", t)
solutionTT = Expression((Hx_code,Hy_code,E_code), t=T, degree=degree)
solutionT = interpolate(solutionTT, All_space)
HT, ET = solutionT.split(deepcopy=True)
# vtkfile1 << ET
# vtkfile3 << HT
Eerror_L2 = errornorm(ET, Esolution, 'L2')#degree-(r-1))
print('Eerror_L2 =',Eerror_L2)
Herror_L2 = errornorm(HT, Hsolution, 'L2')#degree-r)
print('Herror_L2 =',Herror_L2)
Herror_Hcurl = errornorm(HT, Hsolution, 'Hcurl')#degree-r)
print('Herror_Hcurl =',Herror_Hcurl)
errors = {'Eerror_L2' : Eerror_L2, 'Herror_L2' :Herror_L2, 'Herror_Hcurl':Herror_Hcurl}
Error.append(errors)
# Print convergence rates
from math import log as ln
error_types = Error[0].keys()
for error_type in sorted(error_types):
j = len(Error)
print ('\nError norm of', error_type)
for i in range(0, j-1):
r = ln(Error[i+1][error_type]/Error[i][error_type])/ln(0.5) # E is a list of errors
print ('mesh size =%.7f Ei=%.7f Ei+1=%.7f r=%.2f' % (harray[i], Error[i][error_type], Error[i+1][error_type], r))