Code for DG time stepping scheme for heat equation

General
1

Hello everyone,
I am trying to write the code for heat equation using DG method in time and CG in space, applying the algorithm from the book entitled ’ Galerkin Finite Element Methods for Parabolic Problems’ on page no 206, https://www.math.hkust.edu.hk/~mamu/courses/531/parabolic.pdf . Below I am attaching minimal code. Can anyone help me in this regard?

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
#import math
set_log_level(LogLevel.ERROR)
T = 1.0 # final time

# Use SymPy to compute f from the manufactured solution u
import sympy as sym

t = sym.Symbol('t')
x, y = sym.symbols('x[0] x[1]')
#u = sym.exp(-2*t) * sym.sin(x) * sym.sin(y)
u = t*x*(1-x)*y*(1-y)
g = x*y*(1-x)*(1-y)
f = sym.diff(u,t)- sym.diff(sym.diff(u, x), x) - sym.diff(sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)   
print('u =', u_code)
print('f =', f_code) 

for j in range(1,5):
    nx = 2**(j+1)
    T = 1
    num_steps = nx*nx
    dt = 4*T / (num_steps)
    mesh = UnitSquareMesh(nx, nx)
    P1 = FiniteElement('P', triangle, 1)
    element = MixedElement([P1, P1])
    W = FunctionSpace(mesh, element)
    V = FunctionSpace(mesh, 'CG', 1)
    u_D = Expression('t * x[0] * x[1] * (-x[0] + 1) * (-x[1] + 1)', degree=1, t=0)
    u_n = interpolate(u_D, V)
    def boundary(x, on_boundary):
        return on_boundary
    bc1 = DirichletBC(W.sub(0), u_D, boundary)
    bc2 = DirichletBC(W.sub(1), u_D, boundary)
    bc = [bc1, bc2]
    U_0, U_1 = TrialFunctions(W)
    #U_0, U_1 = split(U)
    p, q = TestFunctions(W)
    e_time = 0
    for i in range(num_steps):
        t = (i+1)*dt
        U = Function(W)
        u_D.t = t
        f4 = Expression(f_code, degree=5, t=(i+1)*dt) # final
        f0 = Expression(f_code, degree=5, t=i*dt)
        f1 = Expression(f_code, degree=5, t=(i+ 0.25)*dt)     # initial
        f2 = Expression(f_code, degree=5, t=(i+0.50)*dt)
        f3 = Expression(f_code, degree=5, t=(i+0.75)*dt)
        f = dt/ 12 *( f0 + 4*(f1+f3) + 2*f2 + f4) # using simpson's 1/3rd rulle 
        g0 = -dt*f0 # for right hand side of second equation
        g1 = -0.75*f1
        g2 = -0.50*f2
        g3 = -0.25*f3
        g4 = 0
        g = dt/ 12 *( g0 + 4*(g1+g3) + 2*g2 + g4) # using simpson's 1/3rd rulle

        F = (U_0 + U_1) * p * dx - dt * dot(grad(U_0), grad(p)) * dx - dt / 2 * dot(grad(U_1), grad(p)) * dx + 1 / 2 * U_1 * q * dx - 1 / 2 * dt * dot(grad(U_0), grad(q)) * dx - dt / 3 * dot(grad(U_1), grad(q)) * dx - u_n * p * dx - f * p * dx - 1 / dt * g * q * dx

        
        solve(F ==0, U, bc)
        U_0, U_1 = U.split()
2

Hi Raksha I want to learn dg method can you help me to solve 1d advection equation using fenicsx.

Please share any code if you have any