# 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 \(\Gamma\) is the domain boundary and \(\nu\) is a constant scalar viscosity. The solution \(u\) is sought in some suitable vector-valued function space \(V\). We take the inner product with an arbitrary test function \(v\in V\) and integrate the viscosity term by parts:

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\) is set to a (fairly arbitrary) small constant value:

```
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 \(v\) with
\((u\cdot\nabla)u\), which is the derivative of the velocity in the
direction \(u\). This directional derivative can be written as
`dot(u,nabla_grad(u))`

since `nabla_grad(u)[i,j]`

\(=\partial_i u_j\).
Note once again that for a nonlinear problem, there are no trial functions in
the formulation. These will be created automatically when the residual
is differentiated by the nonlinear solver:

```
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 = File("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 `File`

object manage the decimation. It
supports both interpolation to linears (the default) or projection (by
passing `project_output=True`

when creating the `File`

).
Outputting data is carried out using the `write()`

method
of `File`

objects:

```
outfile.write(project(u, V_out, name="Velocity"))
```

Finally, we loop over the timesteps solving the equation each time and outputting each result:

```
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.