The goal of this chapter is to demonstrate how a range of important PDEs from science and engineering can be quickly solved with a few lines of FEniCS code. We start with the heat equation and continue with a nonlinear Poisson equation, the equations for linear elasticity, the Navier–Stokes equations, and finally look at how to solve systems of nonlinear advection–diffusion–reaction equations. These problems illustrate how to solve time-dependent problems, nonlinear problems, vector-valued problems, and systems of PDEs. For each problem, we derive the variational formulation and express the problem in Python in a way that closely resembles the mathematics.

As a first extension of the Poisson problem from the previous chapter, we consider the time-dependent heat equation, or the time-dependent diffusion equation. This is the natural extension of the Poisson equation describing the stationary distribution of heat in a body to a time-dependent problem.

We will see that by discretizing time into small time intervals and applying standard time-stepping methods, we can solve the heat equation by solving a sequence of variational problems, much like the one we encountered for the Poisson equation.

Our model problem for time-dependent PDEs reads $$ \begin{alignat}{2} {\partial u\over\partial t} &= \nabla^2 u + f \quad &&\hbox{in }\Omega\times(0, T], \tag{3.1}\\ u &= \ub &&\hbox{on } \partial \Omega\times(0, T], \tag{3.2}\\ u &= \uI &&\mbox{at } t=0\tp \tag{3.3} \end{alignat} $$ Here, \( u \) varies with space and time, e.g., \( u=u(x,y,t) \) if the spatial domain \( \Omega \) is two-dimensional. The source function \( f \) and the boundary values \( \ub \) may also vary with space and time. The initial condition \( \uI \) is a function of space only.

A straightforward approach to solving time-dependent PDEs by the finite element method is to first discretize the time derivative by a finite difference approximation, which yields a sequence of stationary problems, and then turn each stationary problem into a variational formulation.

Let superscript \( n \) denote a quantity at time \( t_n \), where \( n \) is an
integer counting time levels. For example, \( u^n \) means \( u \) at time
level \( n \). A finite difference discretization in time first consists
of sampling the PDE at some time level, say \( t_{n+1} \):
$$
\begin{equation}
\left({\partial u \over\partial t}\right)^{n+1} = \nabla^2 u^{n+1} + f^{n+1}\tp
\tag{3.4}
\end{equation}
$$
The time-derivative can be approximated by a difference quotient.
For simplicity and stability reasons, we choose a
simple backward difference:
$$
\begin{equation}
\left({\partial u\over\partial t}\right)^{n+1}\approx {{u^{n+1} - u^n}\over{\dt}},
\tag{3.5}
\end{equation}
$$
where \( \dt \) is the time discretization parameter.
Inserting (3.5) in (3.4) yields
$$
\begin{equation}
{{u^{n+1} - u^n}\over{\dt}} = \nabla^2 u^{n+1} + f^{n+1}\tp
\tag{3.6}
\end{equation}
$$
This is our time-discrete version of the heat equation
(3.1), a so-called *backward Euler* or *implicit
Euler* discretization.

We may reorder (3.6) so that the left-hand side contains the terms with the unknown \( u^{n+1} \) and the right-hand side contains computed terms only. The result is a sequence of spatial (stationary) problems for \( u^{n+1} \), assuming \( u^n \) is known from the previous time step: $$ \begin{align} u^0 &= \uI, \tag{3.7}\\ u^{n+1} - {\dt}\nabla^2 u^{n+1} &= u^n + {\dt} f^{n+1},\quad n=0,1,2,\ldots \tag{3.8} \end{align} $$ Given \( \uI \), we can solve for \( u^0 \), \( u^1 \), \( u^2 \), and so on.

An alternative to (3.8), which can be convenient in implementations, is to collect all terms on one side of the equality sign: $$ \begin{equation} u^{n+1} - {\dt}\nabla^2 u^{n+1} - u^{n} - {\dt} f^{n+1} = 0,\quad n=0,1,2,\ldots \tag{3.9} \end{equation} $$

We use a finite element method to solve (3.7) and either of the equations (3.8) or (3.9). This requires turning the equations into weak forms. As usual, we multiply by a test function \( v\in \hat V \) and integrate second-derivatives by parts. Introducing the symbol \( u \) for \( u^{n+1} \) (which is natural in the program), the resulting weak form arising from formulation (3.8) can be conveniently written in the standard notation: $$ a(u,v)=L_{n+1}(v),$$ where $$ \begin{align} a(u,v) &= \int_\Omega\left(uv + {\dt} \nabla u\cdot \nabla v\right) \dx, \tag{3.10}\\ L_{n+1}(v) &= \int_\Omega \left(u^n + {\dt} f^{n+1}\right)v \dx\tp \tag{3.11} \end{align} $$ The alternative form (3.9) has an abstract formulation $$ F_{n+1}(u;v) = 0,$$ where $$ \begin{equation} F_{n+1}(u; v) = \int_\Omega \left(uv + {\dt} \nabla u\cdot \nabla v - (u^n + {\dt} f^{n+1})v\right) \dx\tp \tag{3.12} \end{equation} $$

In addition to the variational problem to be solved in each time step,
we also need to approximate the initial condition
(3.7). This equation can also be turned into a
variational problem:
$$ a_0(u,v)=L_0(v),$$
with
$$
\begin{align}
a_0(u,v) &= \int_\Omega uv \dx, \tag{3.13}\\
L_0(v) &= \int_\Omega \uI v \dx\tp \tag{3.14}
\end{align}
$$
When solving this variational problem, \( u^0 \) becomes the \( L^2 \)
projection of the given initial value \( \uI \) into the finite element
space. The alternative is to construct \( u^0 \) by just interpolating the
initial value \( \uI \); that is, if \( u^0=\sum_{j=1}^N U^0_j\phi_j \), we
simply set \( U_j=\uI(x_j,y_j) \), where \( (x_j,y_j) \) are the coordinates
of node number \( j \). We refer to these two strategies as computing the
initial condition by either projection or interpolation. Both
operations are easy to compute in FEniCS through a single statement,
using either the `project`

or `interpolate`

function. The most common
choice is `project`

, which computes an approximation to \( \uI \), but in
some applications where we want to verify the code by reproducing
exact solutions, one must use `interpolate`

(and we use such a test
problem here!).

In summary, we thus need to solve the following sequence of variational problems to compute the finite element solution to the heat equation: find \( u^0\in V \) such that \( a_0(u^0,v)=L_0(v) \) holds for all \( v\in\hat V \), and then find \( u^{n+1}\in V \) such that \( a(u^{n+1},v)=L_{n+1}(v) \) for all \( v\in\hat V \), or alternatively, \( F_{n+1}(u^{n+1},v)=0 \) for all \( v\in\hat V \), for \( n=0,1,2,\ldots \).

Our program needs to implement the time-stepping manually, but can rely on FEniCS to easily compute \( a_0 \), \( L_0 \), \( a \), and \( L \) (or \( F_{n+1} \)), and solve the linear systems for the unknowns.

Just as for the Poisson problem from the previous chapter, we construct a test problem that makes it easy to determine if the calculations are correct. Since we know that our first-order time-stepping scheme is exact for linear functions, we create a test problem which has a linear variation in time. We combine this with a quadratic variation in space. We thus take $$ \begin{equation} u = 1 + x^2 + \alpha y^2 + \beta t, \tag{3.15} \end{equation} $$ which yields a function whose computed values at the nodes will be exact, regardless of the size of the elements and \( \dt \), as long as the mesh is uniformly partitioned. By inserting (3.15) into the heat equation (3.1), we find that the right-hand side \( f \) must be given by \( f(x,y,t)=\beta - 2 - 2\alpha \). The boundary value is \( \ub(x, y, t) = 1 + x^2 + \alpha y^2 + \beta t \) and the initial value is \( \uI(x, y) = 1 + x^2 + \alpha y^2 \).

A new programming issue is how to deal with functions that vary in
space *and time*, such as the boundary condition \( \ub(x, y,
t) = 1 + x^2 + \alpha y^2 + \beta t \). A natural solution is to use a
FEniCS `Expression`

with time \( t \) as a parameter, in addition to the
parameters \( \alpha \) and \( \beta \):

```
alpha = 3; beta = 1.2
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
degree=2, alpha=alpha, beta=beta, t=0)
```

This `Expression`

uses the components of `x`

as independent
variables, while `alpha`

, `beta`

, and `t`

are parameters. The
time `t`

can later be updated by

```
u_D.t = t
```

The essential boundary conditions, along the entire boundary in this case, are implemented in the same way as we have previously implemented the boundary conditions for the Poisson problem:

```
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
```

We shall use the variable `u`

for the unknown \( u^{n+1} \) at the new
time step and the variable `u_n`

for \( u^n \) at the previous time
step. The initial value of `u_n`

can be computed by either projection
or interpolation of \( \uI \). Since we set `t = 0`

for the boundary value
`u_D`

, we can use `u_D`

to specify the initial condition:

```
u_n = project(u_D, V)
# or
u_n = interpolate(u_D, V)
```

We may either define \( a \) or \( L \) according to the formulas above, or we may just define \( F \) and ask FEniCS to figure out which terms should go into the bilinear form \( a \) and which should go into the linear form \( L \). The latter is convenient, especially in more complicated problems, so we illustrate that construction of \( a \) and \( L \):

```
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
```

Finally, we perform the time-stepping in a loop:

```
u = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
u_D.t = t
# Solve variational problem
solve(a == L, u, bc)
# Update previous solution
u_n.assign(u)
```

In the last step of the time-stepping loop, we assign the values of
the variable `u`

(the new computed solution) to the variable `u_n`

containing the values at the previous time step. This must be done
using the `assign`

member function. If we instead try to do `u_n = u`

,
we will set the `u_n`

variable to be the *same* variable as `u`

which is not what we want. (We need two variables, one for the values
at the previous time step and one for the values at the current time
step.)

`u_D.t`

must be updated before the
`solve`

statement to enforce computation of Dirichlet conditions at
the current time step. A Dirichlet condition defined in terms of an
`Expression`

looks up and applies the value of a parameter such as `t`

when it gets evaluated and applied to the linear system.
The time-stepping loop above does not contain any comparison of the
numerical and the exact solutions, which we must include in order to
verify the implementation. As for the Poisson equation in
the section Dissection of the program, we compute the difference
between the array of nodal values for `u`

and the array of nodal
values for the interpolated exact solution. This may be done as
follows:

```
u_e = interpolate(u_D, V)
error = np.abs(u_e.vector().array() - u.vector().array()).max()
print('t = %.2f: error = %.3g' % (t, error))
```

For the Poisson example, we used the function
`compute_vertex_values`

to extract the function values at the
vertices. Here we illustrate an alternative method to extract the
vertex values, by calling the function `vector`

, which returns
the vector of degrees of freedom. For a \( \mathsf{P}_1 \)
function space, this vector of degrees of freedom will be equal to
the array of vertex values obtained by calling
`compute_vertex_values`

, albeit possibly in a different order.

The complete program for solving the heat equation goes as follows:

```
from fenics import *
import numpy as np
T = 2.0 # final time
num_steps = 10 # number of time steps
dt = T / num_steps # time step size
alpha = 3 # parameter alpha
beta = 1.2 # parameter beta
# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
degree=2, alpha=alpha, beta=beta, t=0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
u_D.t = t
# Compute solution
solve(a == L, u, bc)
# Plot solution
plot(u)
# Compute error at vertices
u_e = interpolate(u_D, V)
error = np.abs(u_e.vector().array() - u.vector().array()).max()
print('t = %.2f: error = %.3g' % (t, error))
# Update previous solution
u_n.assign(u)
# Hold plot
interactive()
```

This example program can be found in the file `ft03_heat.py`.

Let us now solve a more interesting test problem, namely the diffusion of a Gaussian hill. We take the initial value to be $$ \uI(x,y)= e^{-ax^2 - ay^2}$$ for \( a = 5 \) on the domain \( [-2,2]\times [2,2] \). For this problem we will use homogeneous Dirichlet boundary conditions (\( \ub = 0 \)).

Which are the required changes to our previous program? One major
change is that the domain is no longer a unit square. The new domain can
be created easily in FEniCS using `RectangleMesh`

:

```
nx = ny = 30
mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
```

Note that we have used a much higher resolution than before to better
resolve the features of the solution. We also need to redefine the
initial condition and the boundary condition. Both are easily changed by
defining a new `Expression`

and by setting \( u = 0 \) on the boundary.

To be able to visualize the solution in an external program such as
ParaView, we will save the solution to a file in VTK format in each time
step. We do this by first creating a `File`

with the suffix `.pvd`

:

```
vtkfile = File('heat_gaussian/solution.pvd')
```

Inside the time loop, we may then append the solution values to this file:

```
vtkfile << (u, t)
```

This line is called in each time step, resulting in the creation of
a new file with suffix `.vtu`

containing all data for the time
step (the mesh and the vertex values). The file
`heat_gaussian/solution.pvd`

will contain the time values and
references to the `.vtu`

file, which means that the `.pvd`

file will be a
single small file that points to a large number of `.vtu`

files
containing the actual data. Note that we choose to store the solution
to a subdirectory named `heat_gaussian`

. This is to avoid cluttering
our source directory with all the generated data files.
One does not need to create the directory before running the
program as it will be created automatically by FEniCS.

The complete program appears below.

```
from fenics import *
import time
T = 2.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
# Create mesh and define function space
nx = ny = 30
mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant(0), boundary)
# Define initial value
u_0 = Expression('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))',
degree=2, a=5)
u_n = interpolate(u_0, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
# Create VTK file for saving solution
vtkfile = File('heat_gaussian/solution.pvd')
# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u, bc)
# Save to file and plot solution
vtkfile << (u, t)
plot(u)
# Update previous solution
u_n.assign(u)
# Hold plot
interactive()
```

This example program can be found in the file `ft04_heat_gaussian.py`.

To visualize the diffusion of the Gaussian hill, start ParaView,
choose **File–Open...**, open `heat_gaussian/solution.pvd`

, and click
**Apply** in the Properties pane. Click on the play button to display
an animation of the solution. To save the animation to a file, click
**File–Save Animation...** and save the file to a desired file format,
for example AVI or Ogg/Theora.
Once the animation has been saved to a file, you can play the animation
offline using a player such as mplayer or VLC, or upload your
animation to YouTube. Figure 5 shows a sequence
of snapshots of the solution.