A time-dependent problem, Burgers' equation

We will solve the viscous Burgers equation, a nonlinear equation for the advection and diffusion on momentum in one dimension.

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - \nu \frac{\partial^2 u}{\partial x^2} = 0 $$

We will solve on a periodic interval mesh, and therefore do not impose any boundary conditions. As usual, we need to derive a variational form.

Spatial discretisation

We first discretise in space, mulitplying by a test function $v \in V$ and integrating the viscosity term by parts to obtain the semi-discrete problem. Find $u(x, t) \in V$ such that

$$ \int_\Omega \frac{\partial u}{\partial t} v + u \frac{\partial u}{\partial x} v + \nu \frac{\partial u}{\partial x}\frac{\partial v}{\partial x} \, \mathrm{d}x = 0 \quad \forall v \in V. $$

Time discretisation

We now need to discretise in time. For simplicity, and stability we'll use backward Euler, replacing all instances of $u$ with $u^{n+1}$ and the time derivative by $\frac{u^{n+1} - u^n}{\Delta t}$. We end up with the discrete problem, find $u^{n+1} \in V$ such that

$$ \int_\Omega \frac{u^{n+1} - u^n}{\Delta t} v + u^{n+1} \frac{\partial u^{n+1}}{\partial x} v + \nu \frac{\partial u^{n+1}}{\partial x}\frac{\partial v}{\partial x} \, \mathrm{d}x = 0 \quad \forall v \in V. $$


To solve the problem in a concrete setting, we need two things. A domain, and an initial condition for $u$. For the former, we'll choose a periodic interval of length 2, for the latter, we'll start with $u = \sin(2 \pi x)$.

In addition we need to choose the viscosity, which we will set to a small constant value $\nu = 10^{-2}$

As ever, we begin by importing Firedrake

In [1]:
%matplotlib inline
from firedrake import *

n = 100
mesh = PeriodicIntervalMesh(n, length=2)

x = SpatialCoordinate(mesh)[0]

u_init = sin(2*pi*x)
In [2]:
nu = Constant(1e-2)

We choose degree 2 piecewise continuous Lagrange polynomials for our solution and test space:

In [3]:
V = FunctionSpace(mesh, "Lagrange", 2)

We also need solution functions for $u^{n+1}$ and $u^n$, along with a test function $v$.

In [4]:
u_n1 = Function(V, name="u^{n+1}")
u_n = Function(V, name="u^{n}")
v = TestFunction(V)

We provide the initial condition for $u_n$, and choose a $\Delta t$ such that the advective Courant number is around 1. This is more restrictive than required for stability of the time integration, but gives us enough accuracy to see the temporal evolution of the system.

In [5]:
dt = 1.0 / n

Now we're ready to define the variational form. Since this problem is nonlinear, note that we do not have a trial function anywhere. We just write down the residual, Firedrake will automatically compute the Jacobian by differentiating the residual inside the nonlinear solver.

In [6]:
F = (((u_n1 - u_n)/dt) * v +
     u_n1 * u_n1.dx(0) * v + 

For visualisation purposes, we will save a copy of the state $u_n$ at each timestep, we can plot and animate these in the notebook if the ipywidgets package is installed.

In [7]:
# If passed an existing Function object, the Function 
# constructor makes a copy.
results = [Function(u_n)]

Finally, we loop over the timesteps, solving the equation and advancing in time.

In [8]:
t = 0.0
t_end = 0.5
while t <= t_end:
    solve(F == 0, u_n1)
    t += dt

To visualize the results, we'll create a movie using matplotlib's animation tools. First, we'll create a figure and axes to draw on and plot the initial values.

In [9]:
import matplotlib.pyplot as plt
fig, axes = plt.subplots()
axes.set_ylim((-1., 1.))
plot(results[0], axes=axes)
<matplotlib.patches.PathPatch at 0x7fed08f1cf50>

Next, we'll create a function that tells matplotlib how to draw each frame of the animation, which in our case will just be plotting the value at that timestep. The FuncAnimation function will call this on the list of results that we pass in, together with a given interval in milliseconds between each frame. Finally, we'll use the IPython API to render the animation in the notebook.

In [10]:
from matplotlib.animation import FuncAnimation

def animate(u):
    firedrake.plot(u, axes=axes)
    axes.set_ylim((-1., 1.))

interval = 4e3 * float(dt)
animation = FuncAnimation(fig, animate, frames=results, interval=interval)

from IPython.display import HTML