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