Discontinuous Galerkin Method for Burgers Equation

variational formulation
1

Hello everyone,

I am quite new to Fenics and I’m trying to implement a discontinuous Galerkin scheme to solve some simple PDEs in 1D. I started with the advection equation by mimicking this post Lax-Fridrichs Flux for advection equation which worked well.

Now I am trying to solve Burgers equation both with the Upwind and the Lax-Friedrichs numerical fluxes and I run into some problems. First of all, in the case of the upwind flux, I tried adapting the flux in the link above to the flux in the Burgers equation (basically substituting the advection velocity by the value of the field u which is the jacobian of the flux):

    un = 0.5*(u*n[0] + abs(u*n[0])) 

    flux = jump(un*u)

This approach works for DG elements of order 2 but fails in other cases. Apart from that, I tried writing the upwind flux with a conditional expression which works fine for all cases:

    n0 = n[0]

    flux = conditional(u('-')*n0('+')>0, F(u('+')), F(u('-')))

Any ideas why the first option does not yield the expected solution?

The second question is how should I proceed to implement the Lax-Friedrich flux. In this sense, I followed the indications given in Hesthaven’s book Nodal DG methods, Jan S. Hesthaven where:

   f* = avg(f(u('+'),u('-'))) + 0.5*C * jump(u('+'),u('-'))

Where C is the maximum of the flux jacobian (u in this case). I had two different ideas for the choice of C as:

  • The average of u between the two sides of the facet.
  • The maximum of the value of u between both sides of the facet.

None of these options yield a valid solution.

I attach the full code below, thanks in advance.

from fenics import *
from math import pi
from ufl import shape, rank
import matplotlib.pyplot as plt

#Choice of the numerical flux 'upwind' o 'laxf'
flux_method = 'upwind'

#Generate the mesh 
xL = -1
xR = 1
K = 50
cell_size = (xR-xL)/K
mesh = IntervalMesh(K, xL, xR)

#Parameters for time marching algorithm
T = 0.5
cell_size = (xR-xL)/K
N = 100
dt = T/N 
dt_ufl = Constant(dt)

#Create the function space of order 'ord'
Ord = 4
V = FunctionSpace(mesh, 'DG', Ord)
v = TestFunction(V)
du_trial = TrialFunction(V)

# Initial conditions
u_0 = Expression('-sin(pi*x[0])', degree = 2)
u = interpolate(u_0, V)

#Boundary conditions
def boundaryL(x, on_boundary):
    return on_boundary and near(x[0], xL, DOLFIN_EPS)

uL = Constant(0.0)

def boundaryR(x, on_boundary):
    return on_boundary and near(x[0], xR, DOLFIN_EPS)

uR = Constant(0.0)

bcR = DirichletBC(V, uR, boundaryR)

bcL = DirichletBC(V, uL, boundaryL)

bcs = [bcL,bcR]

#Definition of the flux

n = FacetNormal(mesh)
nu = Constant(0.0)

def F(u):
    return 0.5*u*u - nu*u.dx(0)

#Definition of the numerical fluxes

if flux_method == 'upwind':
    #un = 0.5*(u*n[0] + abs(u*n[0]))
    #flux = jump(un*u)
    n0 = n[0]
    flux = conditional(u('-')*n0('+')>0, F(u('+')), F(u('-')))
elif flux_method == 'laxf':
    flux = avg(F(u)) + 0.5*avg(u)*jump(u)
    #flux = avg(F(u)) + 0.5*conditional(u('+')>=u('-'),u('-'),u('+'))*jump(u)

#Weak Form

a = du_trial*v*dx

L1 = dt_ufl * ( F(u)*v.dx(0)*dx 
              - conditional(F(u)*n[0] <= 0, 0.5*uL*uL*n[0]*v, 0.5*u*u*n[0]*v)*ds
              - jump(v)*flux*dS )


#Time evolution with a SSPRK

from ufl import replace
u1 = Function(V); u2 = Function(V)
L2 = replace(L1, {u: u1}); L3 = replace(L1, {u: u2})

du = Function(V)

t = 0.0
step = 0

while t < T - 0.5*dt:
    
    solve(a == L1, du, bcL)
    u1.assign(u + du)
    
    solve(a == L2, du, bcL)
    u2.assign(0.75*u + 0.25*(u1 + du))
    
    solve(a == L3, du, bcL)
    u.assign((1/3)*u + (2/3)*(u2 + du))
    
    step += 1
    t += dt
    
    if  step % 20 == 0: 
        plot ( u , title = ( ' burgers ' ) )
        plt.grid ( True )
        plt.show()
1 Like
2

Hi DGarciaLahuerta,

I also don’t see the error in your code at the moment but I can very much recomend the dolfin_dg package. It helped me to understand the DGFEM formulation a lot better than when I started with fenics. In the demos folder there is also an example on the burgers equation. Perhaps by comparing you can find your bug :slight_smile:

2 Likes
3

Thanks for the kind review. Here’s the paper (or preprint) for more detail and examples.