DG advection equation with upwinding¶
We next consider the advection equation
in a domain \(\Omega\), where \(\vec{u}\) is a prescribed vector field, and \(q(\vec{x}, t)\) is an unknown scalar field. The value of \(q\) is known initially:
and the value of \(q\) is known for all time on the subset of the boundary \(\Gamma\) in which \(\vec{u}\) is directed towards the interior of the domain:
where \(\Gamma_\mathrm{inflow}\) is defined appropriately.
We will look for a solution \(q\) in a space of discontinuous functions \(V\). A weak form of the continuous equation in each element \(e\) is
where we explicitly introduce the subscript \(e\) since the test functions \(\phi_e\) are local to each element. Using integration by parts on the second term, we get
where \(\vec{n}_e\) is an outward-pointing unit normal.
Since \(q\) is discontinuous, we have to make a choice about how to define \(q\) on facets when we assemble the equations globally. We will use upwinding: we choose the upstream value of \(q\) on facets, with respect to the velocity field \(\vec{u}\). We note that there are three types of facets that we may encounter:
Interior facets. Here, the value of \(q\) from the upstream side, denoted \(\widetilde{q}\), is used.
Inflow boundary facets, where \(\vec{u}\) points towards the interior. Here, the upstream value is the prescribed boundary value \(q_\mathrm{in}\).
Outflow boundary facets, where \(\vec{u}\) points towards the outside. Here, the upstream value is the interior solution value \(q\).
We must now express our problem in terms of integrals over the entire mesh and over the sets of interior and exterior facets. This is done by summing our earlier expression over all elements \(e\). The cell integrals are easy to handle, since \(\sum_e \int_e \cdot \,\mathrm{d}x = \int_\Omega \cdot \,\mathrm{d}x\). The interior facet integrals are more difficult to express, since each facet in the set of interior facets \(\Gamma_\mathrm{int}\) appears twice in the \(\sum_e \int_{\partial e}\). In other words, contributions arise from both of the neighbouring cells.
In Firedrake, the separate quantities in the two cells neighbouring an interior facet are denoted by + and -. These markings are arbitrary – there is no built-in concept of upwinding, for example – and the user is responsible for providing a form that works in all cases. We will give an example shortly. The exterior facet integrals are easier to handle, since each facet in the set of exterior facets \(\Gamma_\mathrm{ext}\) appears exactly once in \(\sum_e \int_{\partial e}\). The full equations are then
As a timestepping scheme, we use the three-stage strong-stability-preserving Runge-Kutta (SSPRK) scheme from [SO88]: to discretise \(\frac{\partial q}{\partial t} = \mathcal{L}(q)\), we set
In this worked example, we reproduce the classic cosine-bell–cone–slotted-cylinder advection test case of [LeV96]. The domain \(\Omega\) is the unit square \(\Omega = [0,1] \times [0,1]\), and the velocity field corresponds to solid body rotation \(\vec{u} = (0.5 - y, x - 0.5)\). Each side of the domain has a section of inflow and a section of outflow boundary. We therefore perform both the inflow and outflow integrals over the entire boundary, but construct them so that they only contribute in the correct places.
As usual, we start by importing Firedrake. We also import the math library to give us access to the value of pi. We use a 40-by-40 mesh of squares.
from firedrake import *
from firedrake.pyplot import FunctionPlotter, tripcolor
import math
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
mesh = UnitSquareMesh(40, 40, quadrilateral=True)
We set up a function space of discontinuous bilinear elements for \(q\), and a vector-valued continuous function space for our velocity field.
V = FunctionSpace(mesh, "DQ", 1)
W = VectorFunctionSpace(mesh, "CG", 1)
We set up the initial velocity field using a simple analytic expression.
x, y = SpatialCoordinate(mesh)
velocity = as_vector((0.5 - y, x - 0.5))
u = Function(W).interpolate(velocity)
Now, we set up the cosine-bell–cone–slotted-cylinder initial condition. The first four lines declare various parameters relating to the positions of these objects, while the analytic expressions appear in the last three lines.
bell_r0 = 0.15; bell_x0 = 0.25; bell_y0 = 0.5
cone_r0 = 0.15; cone_x0 = 0.5; cone_y0 = 0.25
cyl_r0 = 0.15; cyl_x0 = 0.5; cyl_y0 = 0.75
slot_left = 0.475; slot_right = 0.525; slot_top = 0.85
bell = 0.25*(1+cos(math.pi*min_value(sqrt(pow(x-bell_x0, 2) + pow(y-bell_y0, 2))/bell_r0, 1.0)))
cone = 1.0 - min_value(sqrt(pow(x-cone_x0, 2) + pow(y-cone_y0, 2))/cyl_r0, 1.0)
slot_cyl = conditional(sqrt(pow(x-cyl_x0, 2) + pow(y-cyl_y0, 2)) < cyl_r0,
conditional(And(And(x > slot_left, x < slot_right), y < slot_top),
0.0, 1.0), 0.0)
We then declare the initial condition of \(q\) to be the sum of these fields. Furthermore, we add 1 to this, so that the initial field lies between 1 and 2, rather than between 0 and 1. This ensures that we can’t get away with neglecting the inflow boundary condition. We also save the initial state so that we can check the \(L^2\)-norm error at the end.
q = Function(V).interpolate(1.0 + bell + cone + slot_cyl)
q_init = Function(V).assign(q)
Next we’ll create a list to store the function values at every timestep so that we can make a movie of them later.
qs = []
We will run for time \(2\pi\), a full rotation. We take 600 steps, giving
a timestep close to the CFL limit. We declare an extra variable dtc; for
technical reasons, this means that Firedrake does not have to compile new C code
if the user tries different timesteps. Finally, we define the inflow boundary
condition, \(q_\mathrm{in}\). In general, this would be a Function, but
here we just use a Constant value.
T = 2*math.pi
dt = T/600.0
dtc = Constant(dt)
q_in = Constant(1.0)
Now we declare our variational forms. Solving for \(\Delta q\) at each stage, the explicit timestepping scheme means that the left hand side is just a mass matrix.
dq_trial = TrialFunction(V)
phi = TestFunction(V)
a = phi*dq_trial*dx
The right-hand-side is more interesting. We define n to be the built-in
FacetNormal object; a unit normal vector that can be used in integrals over
exterior and interior facets. We next define un to be an object which is
equal to \(\vec{u}\cdot\vec{n}\) if this is positive, and zero if this is
negative. This will be useful in the upwind terms.
n = FacetNormal(mesh)
un = 0.5*(dot(u, n) + abs(dot(u, n)))
We now define our right-hand-side form L1 as \(\Delta t\) times the
sum of four integrals.
The first integral is a straightforward cell integral of
\(q\nabla\cdot(\phi\vec{u})\). The second integral represents the inflow
boundary condition. We only want this to contribute on the inflow part of the
boundary, where \(\vec{u}\cdot\vec{n} < 0\) (recall that \(\vec{n}\) is
an outward-pointing normal). Where this is true, the condition gives the
desired expression \(\phi q_\mathrm{in}\vec{u}\cdot\vec{n}\), otherwise the
condition gives zero. The third integral operates in a similar way to give
the outflow boundary condition. The last integral represents the integral
\(\widetilde{q}(\phi_+ \vec{u} \cdot \vec{n}_+ + \phi_- \vec{u} \cdot \vec{n}_-)\)
over interior facets. We could again use a conditional in order to represent
the upwind value \(\widetilde{q}\) by the correct choice of \(q_+\) or
\(q_-\), depending on the sign of \(\vec{u}\cdot\vec{n_+}\), say.
Instead, we make use of the quantity un, which is either
\(\vec{u}\cdot\vec{n}\) or zero, in order to avoid writing explicit
conditionals. Although it is not obvious at first sight, the expression given in
code is equivalent to the desired expression, assuming
\(\vec{n}_- = -\vec{n}_+\).
L1 = dtc*(q*div(phi*u)*dx
- conditional(dot(u, n) < 0, phi*dot(u, n)*q_in, 0.0)*ds
- conditional(dot(u, n) > 0, phi*dot(u, n)*q, 0.0)*ds
- (phi('+') - phi('-'))*(un('+')*q('+') - un('-')*q('-'))*dS)
In our Runge-Kutta scheme, the first step uses \(q^n\) to obtain
\(q^{(1)}\). We therefore declare similar forms that use \(q^{(1)}\)
to obtain \(q^{(2)}\), and \(q^{(2)}\) to obtain \(q^{n+1}\). We
make use of UFL’s replace feature to avoid writing out the form repeatedly.
q1 = Function(V); q2 = Function(V)
L2 = replace(L1, {q: q1}); L3 = replace(L1, {q: q2})
We now declare a variable to hold the temporary increments at each stage.
dq = Function(V)
Since we want to perform hundreds of timesteps, ideally we should avoid
reassembling the left-hand-side mass matrix each step, as this does not change.
We therefore make use of the LinearVariationalProblem and
LinearVariationalSolver objects for each of our Runge-Kutta stages. These
cache and reuse the assembled left-hand-side matrix. Since the DG mass matrices
are block-diagonal, we use the ‘preconditioner’ ILU(0) to solve the linear
systems. As a minor technical point, we in fact use an outer block Jacobi
preconditioner. This allows the code to be executed in parallel without any
further changes being necessary.
params = {'ksp_type': 'preonly', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}
prob1 = LinearVariationalProblem(a, L1, dq)
solv1 = LinearVariationalSolver(prob1, solver_parameters=params)
prob2 = LinearVariationalProblem(a, L2, dq)
solv2 = LinearVariationalSolver(prob2, solver_parameters=params)
prob3 = LinearVariationalProblem(a, L3, dq)
solv3 = LinearVariationalSolver(prob3, solver_parameters=params)
We now run the time loop. This consists of three Runge-Kutta stages, and every 20 steps we write out the solution to file and print the current time to the terminal.
t = 0.0
step = 0
output_freq = 20
while t < T - 0.5*dt:
solv1.solve()
q1.assign(q + dq)
solv2.solve()
q2.assign(0.75*q + 0.25*(q1 + dq))
solv3.solve()
q.assign((1.0/3.0)*q + (2.0/3.0)*(q2 + dq))
step += 1
t += dt
if step % output_freq == 0:
qs.append(q.copy(deepcopy=True))
print("t=", t)
To check our solution, we display the normalised \(L^2\) error, by comparing to the initial condition.
L2_err = sqrt(assemble((q - q_init)*(q - q_init)*dx))
L2_init = sqrt(assemble(q_init*q_init*dx))
print(L2_err/L2_init)
Finally, we’ll animate our solution using matplotlib. We’ll need to evaluate the solution at many points in every frame of the animation, so we’ll employ a helper class that pre-computes some relevant data in order to speed up the evaluation.
nsp = 16
fn_plotter = FunctionPlotter(mesh, num_sample_points=nsp)
We first set up a figure and axes and draw the first frame.
fig, axes = plt.subplots()
axes.set_aspect('equal')
colors = tripcolor(q_init, num_sample_points=nsp, vmin=1, vmax=2, axes=axes)
fig.colorbar(colors)
Now we’ll create a function to call in each frame. This function will use the helper object we created before.
def animate(q):
colors.set_array(fn_plotter(q))
The last step is to make the animation and save it to a file.
interval = 1e3 * output_freq * dt
animation = FuncAnimation(fig, animate, frames=qs, interval=interval)
try:
animation.save("DG_advection.mp4", writer="ffmpeg")
except:
print("Failed to write movie! Try installing `ffmpeg`.")
This demo can be found as a script in DG_advection.py.
References
Randall J. LeVeque. High-Resolution Conservative Algorithms for Advection in Incompressible Flow. SIAM Journal on Numerical Analysis, 33(2):627–665, 1996. doi:10.1137/0733033.
Chi-Wang Shu and Stanley Osher. Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes. Journal of Computational Physics, 77(2):439–471, 1988. doi:10.1016/0021-9991(88)90177-5.