DG convection convergence rate

Question

Hello!

I'm solving the convection-diffusion equation (without diffusion),
$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = f(x,t) $$
with a 4th order Runge-Kutta scheme in time. I'm testing the code with an MMS-test with exact solution
$$u = \sin (2 \pi x) \sin(2 \pi t) $$
so the source term becomes
$$f(x,t) = 2 \pi \cos (2 \pi t) \sin (2 \pi x) + 2 \pi c \cos (2 \pi x) \sin (2 \pi t). $$

If we eliminate the error in space, I should obtain 4th order convergence in time. However, I only get convergence rate 1, and can't understand why.

Does anyone else understand why?

Thank you!

Regards,
Per Thomas

from dolfin import *
import numpy as np
set_log_level(30)

elements_order = 4
N = 40
T = 0.1

# Initial condition
u0 = Constant(0.0)

# Manufactured solution
exact_ = Expression("sin(2*pi*x[0])*sin(2*pi*t)", t=0.0, degree=elements_order + 3)

# Source term
f = Expression("2*pi*cos(2*pi*t)*sin(2*pi*x[0]) + 2*pi*c*cos(2*pi*x[0])*sin(2*pi*t)", t=0.0, c=0.1, degree=elements_order + 3)

# Velocity
c = Constant(0.1)

def L(u,v,c,mesh,f):
    n = FacetNormal(mesh)
    c_ = as_vector((c,))
    cn = (dot(c_, n) + abs(dot(c_, n)))/2.0
    return inner(c*u,v.dx(0))*dx - dot(cn('+')*u('+') - cn('-')*u('-'), jump(v))*dS - dot(v, cn*u)*ds + f*v*dx

mesh = IntervalMesh(N, 0, 1)

V_DG = FunctionSpace(mesh, "DG", elements_order)
V_DG_higher_order = FunctionSpace(mesh, "DG", elements_order + 3)

u = TrialFunction(V_DG)
v = TestFunction(V_DG)

u_p = Function(V_DG)
u_e = Function(V_DG)
f_ = Function(V_DG)

u_0 = Function(V_DG)
u_1 = Function(V_DG)
u_2 = Function(V_DG)
u_3 = Function(V_DG)
u_4 = Function(V_DG)


err_values = []
dt_values = [0.01, 0.005, 0.001]
for dt in dt_values:

    u_p.assign(u0)
    u_0.assign(u_p)

    F0 = inner(u,v)*dx - u_0*v*dx - dt*L(u_0, v, c, mesh, f)
    F1 = inner(u,v)*dx - u_1*v*dx - dt*L(u_1, v, c, mesh, f)
    F2 = inner(u,v)*dx - u_2*v*dx - dt*L(u_2, v, c, mesh, f)
    F3 = inner(u,v)*dx - 3./8*u_0*v*dx - 1./3*u_1*v*dx - 1./4*u_2*v*dx - 1./24*(u_3*v*dx + dt*L(u_3, v, c, mesh, f))

    t = 0.0
    T = 0.1
    step = 0
    while t < T:
        c.t = t
        f.t = t
        f_.assign(interpolate(f, V_DG))

        u_0.assign(u_p)
        solve(lhs(F0) == rhs(F0), u_1, [])
        solve(lhs(F1) == rhs(F1), u_2, [])
        solve(lhs(F2) == rhs(F2), u_3, [])
        solve(lhs(F3) == rhs(F3), u_4, [])
        u_p.assign(u_4)

        plot(u_p)

        step += 1
        t = step*dt
        print "t =",t

    exact_.t = T
    u_e.assign(interpolate(exact_, V_DG))

    # Computing error
    u_ = interpolate(u_p, V_DG_higher_order)
    e = abs(u_ - u_e)
    L2 = sqrt(assemble(e*e*dx(mesh)))
    err_values.append(L2)

    plot(u_e, title="Exact")

# Convergence rates
m = len(dt_values)
r_values_L2 = [np.log(err_values[i-1]/err_values[i])/
            np.log(dt_values[i-1]/dt_values[i])
            for i in range(1, m, 1)]
print "Error:", err_values
print "Rates:", r_values_L2

interactive()

# Returns:
"""
Error: [0.004292230562314632, 0.0021784458312416196, 0.0004408729573250635]
Rates: [0.97842833431181697, 0.99265104445845187]
"""