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.

The heat equation

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.

PDE problem

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.

Variational formulation

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 \).

FEniCS implementation

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.

Test problem 1: A known analytical solution

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 \).

FEniCS implementation

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

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.

Test problem 2: Diffusion of a Gaussian function

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 \)).

FEniCS implementation

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.

Visualization in ParaView

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.

• «
• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• ...
• 23
• »