Hi everyone,
I’m currently working on a time-dependent nonlinear PDE in FEniCS and performing a temporal convergence study. I’m using a fixed spatial mesh and comparing my numerical solution against a reference solution computed on a much finer time grid.
However, I’m having trouble obtaining a clean optimal convergence rate. The expected temporal order should be around second order, but the observed rate of convergence (ROC) fluctuates quite a bit, especially for coarser time steps.
Here is my output:
dt L2-error L2-order H1-error H1-order
-----------------------------------------------------------
2.000000e-01 1.228031e-05 0.0000 1.091201e-04 0.0000
1.000000e-01 1.441561e-06 3.0906 2.038898e-05 2.4201
5.000000e-02 2.418780e-07 2.5753 2.151336e-06 3.2445
2.500000e-02 5.629508e-08 2.1032 5.005621e-07 2.1036
1.250000e-02 1.338052e-08 2.0729 1.189736e-07 2.0729
-----------------------------------------------------------
The rates eventually settle close to 2, but for larger time steps they fluctuate significantly and even exceed the expected order (e.g., 3.09 and 3.24).
I’m wondering:
- Is this kind of fluctuation normal in temporal convergence studies?
- Could the issue come from the reference solution not being sufficiently fine?
- Is there a possibility that spatial errors are still influencing the temporal error even though I fixed a relatively fine mesh?
I’m attaching my FEniCS code as well for context.
Any suggestions or insights would be greatly appreciated. Thanks in advance!
import numpy as np
from fenics import *
import time
pi = np.pi
T = 10.0
nref = 5
eL2 = []
eH1 = []
hh = []
rL2 = [0.0]
rH1 = [0.0]
solutions = []
spaces = []
meshes = []
total_start = time.time()
def main(N, M, k1):
level_start = time.time()
h = 1.0 / N
dt = 1.0/M
num_steps = int(round(T / dt))
# ── Mesh & space ──
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", k1)
u_0 = Expression(
'sin(2*pi*x[0])*sin(2*pi*x[1])',
degree=2, pi=pi
)
f = Expression(
'-(1-k/2.0)*sin(2*pi*x[0])*sin(2*pi*x[1])',
degree=2,
pi=pi,
k=dt
)
# ── Dirichlet BC ──
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant(0.0), boundary)
# ── Interpolate u^0 ──
u_k = interpolate(u_0, V)
u0 = Function(V)
du0 = TrialFunction(V)
v0 = TestFunction(V)
u0.assign(u_k) # Newton initial guess
# ── First-step residual F0(u0; v0) ───
F0 = (
((u0 - u_k) / dt) * v0 * dx
+ (dt / 4.0) * dot(grad(u0 + u_k), grad(v0)) * dx
+ (dt / 8.0) * (u0**3 + u0**2*u_k + u0*u_k**2 + u_k**3) * v0 * dx
- (dt/2.0) * dot(grad(u0), grad(v0)) * dx -f * v0 * dx
)
J0 = derivative(F0, u0, du0)
solve(
F0 == 0, u0, bc, J=J0,
solver_parameters={
"newton_solver":{
"relative_tolerance":1e-6,
"absolute_tolerance":1e-6,
"maximum_iterations":50,
"linear_solver":"lu"
}
}
)
# ── Initialise 3-level state ───
uoldd = Function(V)
uold = Function(V)
u = Function(V)
uoldd.assign(u_k) # u^{n-1} = u^0
uold.assign(u0) # u^n = u^1
phi = TestFunction(V)
du = TrialFunction(V)
t = dt
for _ in range(1, num_steps):
u.assign(uold) # Newton initial guess
F = (
((u - 2*uold + uoldd) / dt**2) * phi * dx
+ 0.5 * dot(grad(u), grad(phi)) * dx
+ 0.5 * dot(grad(uoldd), grad(phi)) * dx
+ dot(grad(u), grad(phi)) * 1./(2*dt) * dx
- dot(grad(uoldd), grad(phi)) * 1./(2*dt) * dx
+ u * phi * 1./(2*dt) * dx
- uoldd * phi * 1./(2*dt) * dx
+ (1.0/4.0) * u**3 * phi * dx
+ (1.0/4.0) * u**2 * uoldd * phi * dx
+ (1.0/4.0) * uoldd**2 * u * phi * dx
+ (1.0/4.0) * uoldd**3 * phi * dx
)
J = derivative(F, u, du)
solve(
F == 0, u, bc, J=J,
solver_parameters={
"newton_solver":{
"relative_tolerance":1e-6,
"absolute_tolerance":1e-6,
"maximum_iterations":50,
"linear_solver":"lu"
}
}
)
uoldd.assign(uold)
uold.assign(u)
t += dt
return u, V, mesh
if __name__ == "__main__":
degree = 2
Mlist = [5,10,20,40,80]
dtvals = []
eU = []
eUH1 = []
rU = [0.0]
rUH1 = [0.0]
Nref = 100
Mref = 512
print("Computing reference solution")
uh_ref, V_ref, mesh_ref = main(
Nref,
Mref,
degree
)
Nfixed = 100
for M in Mlist:
print("Solving dt = ",1./M)
uh,V,mesh = main(
Nfixed,
M,
degree
)
dt = 1./M
dtvals.append(dt)
# Common space
W = FunctionSpace(
mesh_ref,
"CG",
degree
)
uhW = interpolate(
uh,
W
)
urefW = interpolate(
uh_ref,
W
)
diff = Function(W)
diff.vector()[:] = \
uhW.vector().get_local() \
- urefW.vector().get_local()
# L2 error
L2_error = sqrt(
assemble(
diff*diff*dx
)
)
# H1 error
H1_error = sqrt(
assemble(
inner(
grad(diff),
grad(diff)
)*dx
)
)
eU.append(L2_error)
eUH1.append(H1_error)
k = len(eU)-1
if k>0:
rU.append(
np.log(
eU[k-1]/eU[k]
)/
np.log(
dtvals[k-1]/dtvals[k]
)
)
rUH1.append(
np.log(
eUH1[k-1]/eUH1[k]
)/
np.log(
dtvals[k-1]/dtvals[k]
)
)
print("\n-----------------------------------------------------------")
print("dt L2-error L2-order H1-error H1-order")
print("-----------------------------------------------------------")
for i in range(len(eU)):
print(
"{:.6e} {:.6e} {:.4f} {:.6e} {:.4f}".format(
dtvals[i],
eU[i],
rU[i],
eUH1[i],
rUH1[i]
)
)
print("-----------------------------------------------------------")