Behaviour of 1D discontinuous Galerkin with LF flux for Burgers' eq

Hello all. I am new to both fenics and finite element problems in general. I am trying to replicate the results found here: Bitbucket using implicit Euler time-stepping rather than spacetime finite elements.

However, my code as constructed causes the appearance of an artefact near the shock that is resistant to changes in the size of the timestep.

I’m unsure if this problem is caused by the change in approach to the time dimension or an incorrect specification of the Lax-Friedrichs flux

import matplotlib.pyplot as plt
import ufl
from dolfin import (
    Constant, IntervalMesh, SubDomain, DOLFIN_EPS, project, inner,
    as_vector, PETScSNESSolver, dot, dS, NonlinearProblem, plot,
    FunctionSpace, TrialFunction, TestFunction, Expression, MeshFunction,
    Measure, FacetNormal, derivative, assemble, dx, interpolate, Function,
    solve, DirichletBC, avg, jump, nabla_grad, grad)
import time
import numpy as np

start = time.time()

#domain parameters
x_R = 1.
x_L = 0.
N_el = 64
L = x_R - x_L
cell_size = (x_R-x_L)/N_el

#Parameters for time marching algorithm
T = 0.5
dt = 0.5/64
dt_ufl = Constant(dt)

nu = Constant(0.01)

# Mesh and function space.
mesh = IntervalMesh(N_el, x_L, x_R) # mesh
V = FunctionSpace(mesh, 'DG', 1) #function space

# Boundary and initial conditions
# Set up IC
u0_expr = Expression('max(sin(2.0*pi*x[0] ), 0.0 )', element=V.ufl_element())
Dirichelet_inlet = Constant(0) # Dirichlet BC at inlet

#functions to define the boundary subdomains
class Inlet_BC(SubDomain):
    def inside(self, x, on):
        return (abs(x[0] - x_L) < DOLFIN_EPS) and on

class Outlet_BC(SubDomain):
    def inside(self, x, on):
        return (abs(x[0] - x_R) < DOLFIN_EPS) and on
class Burgers_Variational_form():
    def __init__(self, x_L, x_R, N_el, dt, nu, mesh, V, u0, Inlet_BC):
        self.x_L = x_L
        self.x_R = x_R
        self.N_el = N_el
        self.dt = dt = nu
        self.dt_ufl = Constant(dt)

        # Mesh and function space.
        self.mesh = mesh
        self.V = V
        self.u = Function(V)
        self.phi = TestFunction(V)
        self.n = FacetNormal(mesh)

        #BCs and ICs
        self.u_old = interpolate(u0, V)
        # self.u.assign(self.u_old)

        # Boundary conditions
        exterior_bdries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
        inlet_ds = Inlet_BC()
        inlet_ds.mark(exterior_bdries, 1)
        outlet_ds = Outlet_BC()
        outlet_ds.mark(exterior_bdries, 2)
        self.ds = Measure('ds', domain=mesh, subdomain_data=exterior_bdries)

    def Fc(self, u):
        flux = u**2/2
        return flux

    # Define local-Lax Friedrichs flux for the convective term
    def numerical_cflux(self, u_p, u_m, n):
        # note that n^+ = -n^-
        C = ufl.Max(abs(u_p), abs(u_m))
        # flux = {{fc(u)}} + C/2*[[u]]
        flux = (self.Fc(u_p) + self.Fc(u_m))/2 + (C/2)*(u_p*n[0] - u_m*n[0])
        return flux

    def domain_numerical_flux(self, u, n):
        interior_numerical_fluxes = inner(self.numerical_cflux(u('+'), u('-'), n('+')), self.phi('+') - self.phi('-'))*dS
        # Exterior Dirichlet condition
        inlet_flux_BC = self.numerical_cflux(u, Dirichelet_inlet, n)*self.phi*self.ds(1)
        # Exterior Neumann condition at outlet
        outlet_flux_BC = self.numerical_cflux(u, u, n)*self.phi*self.ds(2)
        return interior_numerical_fluxes + inlet_flux_BC + outlet_flux_BC

    #volumetric term
    def volume_term(self, u, u_old):
        volume = (inner(u - u_old, self.phi) - self.dt_ufl * inner(self.Fc(u), self.phi.dx(0) )) * dx
        return volume
    def get_residual(self):
        residual = self.volume_term(self.u, self.u_old) + self.dt_ufl * self.domain_numerical_flux(self.u, self.n)
        J = derivative(residual, self.u)
        return residual, J

if __name__ == '__main__':
    form = Burgers_Variational_form(x_L, x_R, N_el, dt, nu, mesh, V, u0_expr, Inlet_BC)
    residual, J = form.get_residual()
    u = form.u
    u_old = form.u_old

    class BurgersProblem(NonlinearProblem):
        def F(self, b, x):
            assemble(residual, tensor=b)
        def J(self, A, x):
            assemble(J, tensor=A)

    burgers = BurgersProblem()
    solver = PETScSNESSolver('newtonls')

    # Time stepping
    t = 0
    i = 0
    while t < T - 0.5*dt:
        solver.solve(BurgersProblem(), u.vector())
        t += dt
        i += 1

A comparison between the example solution and my solution can be seen below:

Any thoughts or advice would be greatly appreciated. Thank you!

At a glance your results look pretty decent. DG methods are still prone to Gibbs phenomenon at shocks. There’s a wealth of literature on flux limiting and shock capturing methods which may help. Also you can see what happens at the shock as you refine the mesh and time step.

I vaguely remember I wrote an example comparing flux approximation techniques for this problem. I’ll see if I can find it…

1 Like

Thank you Nate for your kind response. The response seems to improve with increasing fidelity in the spatial dimension. It worsens with increased temporal resolution, although I have no intuition as to why this might be as it seems counter to the CFL condition.

Thank you also for your work on the Dolfin-dg demos. I’ve found them to be a very useful resource!

If you have any constructive criticism or “wish-list” requests for dolfin_dg, particularly related to DOLFINx, please let me know. I don’t get much time for it anymore, but it could definitely use improvement.