# Discontinuous Galerkin Method for Burgers Equation

**URL:** <https://fenicsproject.discourse.group/t/discontinuous-galerkin-method-for-burgers-equation/5508>\
**Category:** variational formulation\
**Created:** [April 5, 2021, 12:01pm UTC](https://fenicsproject.discourse.group/t/discontinuous-galerkin-method-for-burgers-equation/5508 "2021-04-05T12:01:07Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![DGarciaLahuerta](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dgarcialahuerta/32/2673_2.png) [@DGarciaLahuerta](https://fenicsproject.discourse.group/u/DGarciaLahuerta)\
**Post date:** [April 5, 2021, 12:01pm UTC](https://fenicsproject.discourse.group/t/discontinuous-galerkin-method-for-burgers-equation/5508/1 "2021-04-05T12:01:07Z")

</div>

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](https://fenicsproject.discourse.group/t/lax-friedrichs-flux-for-advection-equation/4647) 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):

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

```

```auto
    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:

```auto
    n0 = n[0]

```

```auto
    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](https://link.springer.com/book/10.1007/978-0-387-72067-8) where:

```auto
   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.

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![RebeccaE](https://avatars.discourse-cdn.com/v4/letter/r/a4c791/32.png) [@RebeccaE](https://fenicsproject.discourse.group/u/RebeccaE)\
**Post date:** [May 6, 2021, 10:26am UTC](https://fenicsproject.discourse.group/t/discontinuous-galerkin-method-for-burgers-equation/5508/2 "2021-05-06T10:26:54Z")

</div>

Hi DGarciaLahuerta,

I also don’t see the error in your code at the moment but I can very much recomend the [dolfin\_dg](https://bitbucket.org/nate-sime/dolfin_dg/src/master/demo/dolfin/burgers/burgers_spacetime.py) 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 🙂

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [May 6, 2021, 2:03pm UTC](https://fenicsproject.discourse.group/t/discontinuous-galerkin-method-for-burgers-equation/5508/3 "2021-05-06T14:03:38Z")

</div>

Thanks for the kind review. Here’s the [paper](https://epubs.siam.org/doi/10.1137/17M1129751) (or [preprint](https://arxiv.org/pdf/1804.02338.pdf)) for more detail and examples.