Burgers equation¶
The Burgers equation is a non-linear equation for the advection and diffusion of momentum. Here we choose to write the Burgers equation in two dimensions to demonstrate the use of vector function spaces:
where
The boundary condition has been used to discard the surface integral. Next, we need to discretise in time. For simplicity and stability we elect to use a backward Euler discretisation:
We can now proceed to set up the problem. We choose a resolution and set up a square mesh:
from firedrake import *
n = 30
mesh = UnitSquareMesh(n, n)
We choose degree 2 continuous Lagrange polynomials. We also need a piecewise linear space for output purposes:
V = VectorFunctionSpace(mesh, "CG", 2)
V_out = VectorFunctionSpace(mesh, "CG", 1)
We also need solution functions for the current and the next timestep. Note that, since this is a nonlinear problem, we don’t define trial functions:
u_ = Function(V, name="Velocity")
u = Function(V, name="VelocityNext")
v = TestFunction(V)
For this problem we need an initial condition:
x = SpatialCoordinate(mesh)
ic = project(as_vector([sin(pi*x[0]), 0]), V)
We start with current value of u set to the initial condition, but we also use the initial condition as our starting guess for the next value of u:
u_.assign(ic)
u.assign(ic)
nu = 0.0001
The timestep is set to produce an advective Courant number of around 1. Since we are employing backward Euler, this is stricter than is required for stability, but ensures good temporal resolution of the system’s evolution:
timestep = 1.0/n
Here we finally get to define the residual of the equation. In the advection
term we need to contract the test function dot(u,nabla_grad(u))
since nabla_grad(u)[i,j]
F = (inner((u - u_)/timestep, v)
+ inner(dot(u,nabla_grad(u)), v) + nu*inner(grad(u), grad(v)))*dx
We now create an object for output visualisation:
outfile = VTKFile("burgers.pvd")
Output only supports visualisation of linear fields (either P1, or
P1DG). In this example we project to a linear space by hand. Another
option is to let the VTKFile
object manage the
decimation. It supports both interpolation to linears (the default) or
projection (by passing project_output=True
when creating the
VTKFile
). Outputting data is carried out using
the write()
method of
VTKFile
objects:
outfile.write(project(u, V_out, name="Velocity"))
Finally, we loop over the timesteps solving the equation each time and outputting each result. Firedrake’s default solver parameters are used, which amount to applying a full LU decomposition as a preconditioner.
t = 0.0
end = 0.5
while (t <= end):
solve(F == 0, u)
u_.assign(u)
t += timestep
outfile.write(project(u, V_out, name="Velocity"))
A python script version of this demo can be found here.