Did I implement this simple method correctly?

[APaganini]

Dear firedrakers,
I tried to implement a simple DG method for the simple advection equation u_t + u_x = 0 in (0,1) with the Dirichlet boundary condition u(t,x=0)=0 and an initial condition u(t=0,x)=u_0(x), but it doesn't work with DG elements (how sad :( )

Could you please confirm whether the firedrake code listed at the end is a correct translation of the following method, which is based on DG+Heun's method:

1. u_ = u(t) + dt*f(u(t))
2. u(t+dt) = u(t) + dt(f(u(t)) + f(u_))/2

where f(u) = -M^{-1}*L(u), with Mdenoting a mass matrix and L(u) = \int u_x v dx + \sum_j [u(xj+)-u(xj-)]v(xj+) (where u(xj+) denotes the value of u at the node xj from the right).

Here is the code. It works with CG elements, but not with DG ones :( (bonus question: DirichletBC(V, Constant(0.), 1) doesn't seem to do anything with DG0 elements, is this correct? I thought DG elements were extended to the bdry,)

from firedrake import *
import matplotlib.pyplot as plt 

# spatial and time discretization
n = 100; dt = 1.0/(10*n) #overkill CFL to be safe

# initialize mesh and set initial condition
mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, "CG", 1) #change to DG to test
u = Function(V) #solution
x = SpatialCoordinate(mesh)[0]
u.interpolate((10**4*(0.1-x)**2*(0.2-x)**2)*(x<0.2)*(x>0.1))

# set up Heun's method
phi, v = TrialFunction(V), TestFunction(V)
a = phi * v * dx
L1 = -(Constant(0.)+u.dx(0))*v*dx - (u('+')-u('-'))*v('+')*dS
u_ = Function(V) #soln at intermediate time
L2 = replace(L1, {u: u_})
fu, fu_ = Function(V), Function(V)
bcs = DirichletBC(V, Constant(0.), 1)
prob1 = LinearVariationalProblem(a, L1, fu, bcs=bcs)
solv1 = LinearVariationalSolver(prob1)
prob2 = LinearVariationalProblem(a, L2, fu_, bcs=bcs)
solv2 = LinearVariationalSolver(prob2)

fig, axes = plt.subplots()
for ii in range(500):
    solv1.solve()
    u_.assign(u + dt*fu)
    solv2.solve()
    u.assign(u + dt/2*(fu+fu_))
    if ii % 100 == 0:
        plot(u, axes=axes)
plt.show()

[dham]

I think you have misunderstood the meaning of the "+" restriction in UFL. "+" and "-" are fixed but arbitrary sides of each facet. They do not have a global geometric meaning.

[APaganini]

Many thanks for pointing this out, David. I was misled by my own notation. I'll fix it and post the corrected version below.

[danshapero]

I've always worked under the assumption that the forms should be the same if you switch all the '+' to '-' and vice versa in dS integrals. Is that correct or have I misunderstood?

[wence-]

I've always worked under the assumption that the forms should be the same if you switch all the '+' to '-' and vice versa in dS integrals. Is that correct or have I misunderstood?

No, that's right, your form should be the same under arbitrary relabelling of facet sides (since a different numbering will produce a different labelling)

[jshipton]

There is an example in the documentation that describes the use of "+" and "-" for the DG advection equation - see here: https://www.firedrakeproject.org/demos/DG_advection.py.html

[APaganini]

Thanks you Jemma, the example there is slightly different, but it's useful to find out what I should do instead.

[APaganini]

As David correctly pointed out, the issue was my my misunderstanding of UFL's "+"and "-". Here is the correct implementation. To experiment with it, you can change the finite element space to "DG0", "DG1", "DG2", "CG1", "CG2", etc.

from firedrake import *
import matplotlib.pyplot as plt

# spatial and time discretization
n = 100; dt = 1.0/(10*n) #overkill CFL to be safe

# initialize mesh and set initial condition
mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, "DG", 1)
u = Function(V) # solution
x = SpatialCoordinate(mesh)[0]
u.interpolate((10**4*(0.1-x)**2*(0.2-x)**2)*(x<0.2)*(x>0.1))

# set up Heun's method
phi, v = TrialFunction(V), TestFunction(V)
a = phi * v * dx
Lv = (Constant(0.)+u.dx(0))*v*dx # volume term
# internal edges term
n = dot(Constant(as_vector([1.])), FacetNormal(mesh))
fr = 0.5*(abs(n) - n) #fr = 1 if n = -1, else fr = 0
Lb = n('-')*(u('+')-u('-'))*(fr('+')*v('+')+fr('-')*v('-'))*dS
L1 = - Lv - Lb
u_ = Function(V) #soln at intermediate time
L2 = replace(L1, {u: u_})
fu, fu_ = Function(V), Function(V)
bcs = DirichletBC(V, Constant(0.), 1)
prob1 = LinearVariationalProblem(a, L1, fu, bcs=bcs)
solv1 = LinearVariationalSolver(prob1)
prob2 = LinearVariationalProblem(a, L2, fu_, bcs=bcs)
solv2 = LinearVariationalSolver(prob2)

fig, axes = plt.subplots()
for ii in range(500):
    solv1.solve()
    u_.assign(u + dt*fu)
    solv2.solve()
    u.assign(u + dt/2*(fu+fu_))
    if ii % 100 == 0:
        plot(u, axes=axes)
plt.show()

[wence-]

bonus question: DirichletBC(V, Constant(0.), 1) doesn't seem to do anything with DG0 elements, is this correct? I thought DG elements were extended to the bdry,

There is no topological association with the boundary facets. we used to support setting things by selecting those basis functions with geometric support on the boundary, but removed it because it's mathematically ill-defined: changing the basis can change which nodes are considered to be "boundary nodes" with a geometric definition.

If you want inflow/outflow conditions on a DG space, enforce them weakly.