Getting order of convergence wrong while applying DG(1)CG(2) method for heat equation

General
1

If I use polynomial degree 2, then while using DG(1)CG(2), i.e., DG time stepping with standard Galerkin method in space, I am not getting the correct order of convergence in L2 and H1 norm. One thing here to mention is that while defining the solution U(t)=~U0+~U1(t−tn−1)/Δt on Jn:=(tn−1,tn], so I think the issue is in defining the initial value because tn−1 is not in Jn, the minimal heat equation code for dg time stepping method for approximating polynomial is given below. Here, as I decrease the mesh size, the error increases. This is minimal code for DG(1) in time; DG(0) reduces to the modified Eulers method in time.

scheme
scheme1533×779 79.2 KB

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

## Equation to solve u_t - (1/2pi^2)\Delta u = 0, for exact solution u_exact = e^{-t} sin(pi*x) sin(pi*y)

set_log_level(LogLevel.ERROR)
T = 1.0 # final time
e_time = 0 
for j in range(1,5):
    nx = 2**(j+1)
    T = 1
    num_steps = nx
    dt = T / (num_steps)
    mesh = UnitSquareMesh(nx, nx)
    P2 = FiniteElement("Lagrange", triangle, 2)
    #element = VectorElement("P2", triangle, 2, dim=2)
    element = MixedElement([P2, P2])
    V = FunctionSpace(mesh, element)
    W = V.sub(0).collapse()
    R = Function(W)
    t = 0
    u_D = Expression('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', degree = 3, t=t)

    def boundary(x, on_boundary):
        return on_boundary
    bc1 = DirichletBC(V.sub(0), u_D, boundary)
    bc2 = DirichletBC(V.sub(1), u_D, boundary)
    bc = [bc1, bc2]
    
    U_1, U_2 = TrialFunctions(V) 
    u_g = Expression(('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', '0'), degree = 3, t=t)
    U_n = interpolate(u_g, V)
    U_n1, U_n2 = split(U_n)
    p, q = TestFunctions(V)
    e_time = 0
    for i in range(num_steps):
        (U_1, U_2) = TrialFunctions(V)
        t = (i+1)*dt
        u_D = Expression('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', degree = 3, t=(i+1)*dt)
        a = (U_1 + U_2) * p * dx + dt/(2*pi*pi) * dot(grad(U_1), grad(p)) * dx + dt / (4*pi*pi) * dot(grad(U_2), grad(p)) * dx + 1 / 2 * U_2 * q * dx + 1 / (4*pi*pi) * dt * dot(grad(U_1), grad(q)) * dx + dt / (6*pi*pi) * dot(grad(U_2), grad(q)) * dx 
        L =  U_n1 * p * dx + U_n2*p*dx 
        U = Function(V)
        solve(a == L, U, bc)
        U_1, U_2 = U.split(deepcopy=True)  
        R = project(U_1 + U_2, W) 
        U_1.assign(U_1)
        U_2.assign(U_2)
    
    vertex_values_u_D = u_D.compute_vertex_values(mesh)
    vertex_values_R = R.compute_vertex_values(mesh)
    error_max = np.max(np.abs(vertex_values_u_D - vertex_values_R))  
    
    error_L2 = errornorm(u_D, R, 'L2')
    error_H1 = errornorm(u_D, R, 'H1')
    e_time = (e_time + dt*error_H1*error_H1)
    e_time1 = sqrt(e_time)
    print ('dt=%.6f nx=%.2f error_L2=%8.4E error_H1=%8.4E e_time1=%8.4E error_max = %8.4E' % (dt, nx, error_L2, error_H1, e_time1, error_max ))

getting wrong results as:

dt=0.250000 nx=4.00 error_L2=2.0440E-01 error_H1=9.3619E-01 e_time1=4.6810E-01 error_max = 4.1447E-01
dt=0.125000 nx=8.00 error_L2=2.5724E-01 error_H1=1.1719E+00 e_time1=4.1431E-01 error_max = 5.1490E-01
dt=0.062500 nx=16.00 error_L2=2.8576E-01 error_H1=1.3014E+00 e_time1=3.2535E-01 error_max = 5.7155E-01
dt=0.031250 nx=32.00 error_L2=3.0068E-01 error_H1=1.3693E+00 e_time1=2.4206E-01 error_max = 6.0136E-01
2

This code is all over the place redefining U_1 and U_2, and never assigning anything to U_n1, U_n2.

3

I am changing it to

       assign(U_n1, U_1)
       assign(U_n2, U_2)

Now getting error Traceback (most recent call last): File "degg_1,2.py", line 47, in <module> assign(U_n1, U_1) AttributeError: 'Indexed' object has no attribute '_cpp_object.'
To avoid assign function if I append like this

    e_time = 0
    UF = [U_n1]
    US = [U_n2]
    for i in range(num_steps):
        (U_1, U_2) = TrialFunctions(V)
        UF.append(U_1)
        US.append(U_2)
        t = (i+1)*dt
        u_D = Expression('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', degree = 3, t=(i+1)*dt)
        F = (U_1 + U_2) * p * dx + dt/(2*pi*pi) * dot(grad(U_1), grad(p)) * dx + dt / (4*pi*pi) * dot(grad(U_2), grad(p)) * dx + 1 / 2 * U_2 * q * dx + 1 / (4*pi*pi) * dt * dot(grad(U_1), grad(q)) * dx + dt / (6*pi*pi) * dot(grad(U_2), grad(q)) * dx -   UF[i] * p * dx - US[i]*p*dx
        
        a, L = lhs(F), rhs(F)
        U = Function(V)
        solve(a == L, U, bc)
        U_1, U_2 = U.split(deepcopy=True)  
        R = project(U_1 + U_2, W)

still getting results as

dt=0.250000 nx=4.00 error_L2=1.8394E-01 error_H1=8.3767E-01 e_time1=4.1883E-01 error_max = 3.6788E-01
dt=0.125000 nx=8.00 error_L2=1.8394E-01 error_H1=8.3767E-01 e_time1=2.9616E-01 error_max = 3.6788E-01
dt=0.062500 nx=16.00 error_L2=1.8394E-01 error_H1=8.3767E-01 e_time1=2.0942E-01 error_max = 3.6788E-01
dt=0.031250 nx=32.00 error_L2=1.8394E-01 error_H1=8.3767E-01 e_time1=1.4808E-01 error_max = 3.6788E-01
4

To check whether you get the theoretical convergence rates, you could also modify my (academic) code that has been created to explain dG(r)cG(s) discretizations applied to a nonlinear heat equation for students. David Kamensky has also some code online for dG(1) time discretizations LinearDGSpaceTimeIntegrator, but I have not tried it out. But basically since you have a linear problem you could also assemble everything as a Kronecker product of temporal and spatial finite element matrices.

2 Likes
5

Thanks, Julian Roth, for sharing the links. I will check it.

1 Like