integration of time-dependent PDE with higher order schemes

Question

Hello,
Following the example of Cahn-Hilliard, I tried to implement for a time-dependent PDE the backward Euler and Trapezoidal rule in time

mesh = IntervalMesh(M, 0, 1)
V = FunctionSpace(mesh, 'CG', 1)
...
# ---------- Euler
z_old = Function(V)
z_old = interpolate(u_start,V)
z = Function(V)
z = interpolate(u_start,V)
dz = TrialFunction(V)
w = TestFunction(V)

F_N = inner(z-z_old,w)*dx + dt*dot(grad(z), grad(w))*dx - inner(dt*rhs(z),w)*dx
dF_N = derivative(F_N, z, dz)

problem = NonlinearVariationalProblem(F_N, z, bcs, J = dF_N)
pdesys_newton = NonlinearVariationalSolver(problem)

# ---------- CN
z2_old = interpolate(u_start,V)
z2 = Function(V)
z2 = interpolate(u_start,V)
dz2 = TrialFunction(V)
w2 = TestFunction(V)

z2_cn = 1/2*z2 + 1/2*z2_old
rhs_cn = 1/2*(rhs(z2) + rhs(z2_old))

F_N2 = z2*w2*dx - z2_old*w2*dx  + dt*dot( grad(z2_cn), grad(w2))*dx - (dt*rhs_cn)*w2*dx
dF_N2 = derivative(F_N2, z2, dz2)

problem_CN = NonlinearVariationalProblem(F_N2, z2, bcs, J = dF_N2)
pdesys_newton_CN = NonlinearVariationalSolver(problem_CN)

but the $L^2$-Error using the Trapezoidal rule is nearly two times the error using the backward Euler in time, and I have no idea why. Could you please give my some suggestions what's going wrong or where to look at? For integration in time I use $\Delta t = \frac{1}{M}$ as step size.

Answer 1

the problem consisted of two parts:

• my rhs was time dependent, so I had to define rhs_old, for which I updated the time step t accordingly
• I used the wrong norm, error = errornorm(u_exact, z2,'l2',5) is correct, as numpy.linalg.norm is a norm for vectors, and not functions.

My CN - method reads now:

# ---------- CN
z2_old = interpolate(u_start,V)
z2 = Function(V)
z2 = interpolate(u_start,V)
dz2 = TrialFunction(V)
w2 = TestFunction(V)

z2_cn = 1/2*z2 + 1/2*z2_old
rhs_cn = 1/2*(rhs(z2) + rhs_old(z2_old))

F_N2 = z2*w2*dx - z2_old*w2*dx  + dt*dot( grad(z2_cn), grad(w2))*dx - (dt*rhs_cn)*w2*dx
dF_N2 = derivative(F_N2, z2, dz2)

problem_CN = NonlinearVariationalProblem(F_N2, z2, bcs, J = dF_N2)
pdesys_newton_CN = NonlinearVariationalSolver(problem_CN)

Answer 2

You should take a look on Gryphon project. It provides different higher-order time integration schemes based on Runge-Kutta familly. Error estimation and time step adaptivity are included too.
https://launchpad.net/gryphonproject

Comments

Yesterday while looking for another possibility I found out about it. But as I'm implementing IMEX-Methods, I have to, at least to a certain extent, implement the methods on my own.