I am trying to implement a DG solver for the linear scalar wave equation from the book of Warburton and Hesthaven (Nodal Discontinuous Galerkin Methods, p20).
So the 1D equation is
$$
\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0
$$
The proposed numerical flux is
$$
\left(au\right)^* = \left\{\left\{au\right\}\right\} + \left|a\right|\frac{1-\alpha}{2}\left[u\right], \qquad 0 \leq \alpha \leq 1
$$
But when I try to implement this, I get a shape mismatch in my numerical flux.
# Define fluxes on interior and exterior facets
uhat = avg(uold) + a*(1.-alpha)/2.0*jump(uold, n)
uhatbnd = uold + a*(1.-alpha)/2.0*uold*n
Does anyone know what's wrong and how to fix it?
My code:
from dolfin import *
# constants
a = 2.0*DOLFIN_PI
lmbda = 0.5*a
l = a/lmbda
omega = -a*a/lmbda
# prepare a mesh
mesh = IntervalMesh(50, 0.0, a)
# Define function spaces
DGorder = 1
V = FunctionSpace(mesh, "DG", DGorder)
U = FunctionSpace(mesh, "DG", DGorder)
# Set initial values
uinit = Expression("0")
u0 = interpolate(uinit, U)
#Set boundary conditions
ubdr = Expression("sin(l*x[0]-omega*t)", l=l, omega=omega, t=0.0)
Port = CompiledSubDomain("x[0] == 0.0")
End = CompiledSubDomain("x[0] == a", a = a)
bdr = [DirichletBC(V, ubdr, Port)]#, DirichletBC(V, uinit, End)]
# Define test and trial functions
v = TestFunction(V)
u = TrialFunction(U)
# Set time stepping parameters
t = 0
dt = 0.01
T = DOLFIN_PI
# Define normals
n = FacetNormal(mesh)
# Set penalty parameter
alpha = 1.
# Central differences for time and spatial discretisation
udot = (u-u0)/dt
uold = (u+u0)*0.5
# Define fluxes on interior and exterior facets
uhat = avg(uold) + a*(1.-alpha)/2.0*jump(uold, n)
uhatbnd = uold + a*(1.-alpha)/2.0*uold*n
# Define variational formulation
F = (udot*v - a*(uold*v.dx(0)))*dx \
- a*inner(uhat, jump(v,n))*dS \
- a*inner(uhatbnd, v*n)*ds
a, L = lhs(F), rhs(F)
# Prepare solution
u = Function(U)
plot(u0)
# the acutal timestepping
while t < T:
solve(a == L, u, bdr)
u0.assign(u)
t += dt
plot(u0)
