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

which yields a function whose computed values at the nodes will be exact, regardless of the size of the elements and \({\Delta t}\), as long as the mesh is uniformly partitioned. By inserting (31) into the heat equation (17), we find that the right-hand side \(f\) must be given by \(f(x,y,t)=\beta - 2 - 2\alpha\). The boundary value is \(u_{_\mathrm{D}}(x, y, t) = 1 + x^2 + \alpha y^2 + \beta t\) and the initial value is \(u_{_0}(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 \(u_{_\mathrm{D}}(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 \(u_{_0}\). 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

for \(a = 5\) on the domain \([-2,2]\times [2,2]\). For this problem we will use homogeneous Dirichlet boundary conditions (\(u_{_\mathrm{D}} = 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 A sequence of snapshots of the solution of the Gaussian hill problem created with ParaView shows a sequence of snapshots of the solution.

A sequence of snapshots of the solution of the Gaussian hill problem created with ParaView

Test problem¶

To solve a test problem, we need to choose the right-hand side \(f\), the coefficient \(q(u)\) and the boundary value \(u_{_\mathrm{D}}\). Previously, we have worked with manufactured solutions that can be reproduced without approximation errors. This is more difficult in nonlinear problems, and the algebra is more tedious. However, we may utilize SymPy for symbolic computing and integrate such computations in the FEniCS solver. This allows us to easily experiment with different manufactured solutions. The forthcoming code with SymPy requires some basic familiarity with this package. In particular, we will use the SymPy functions diff for symbolic differentiation and ccode for C/C++ code generation.

We take \(q(u) = 1 + u^2\) and define a two-dimensional manufactured solution that is linear in \(x\) and \(y\):

# Warning: from fenics import * will import both `sym` and
# `q` from FEniCS. We therefore import FEniCS first and then
# overwrite these objects.
from fenics import *

def q(u):
    "Return nonlinear coefficient"
    return 1 + u**2

# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(q(u)*sym.diff(u, x), x) - sym.diff(q(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

Turning the expressions for u and f into C or C++ syntax for FEniCS Expression objects needs two steps. First, we ask for the C code of the expressions:

u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)

In some cases, one will need to edit the result to match the required syntax of Expression objects, but not in this case. (The primary example is that M_PI for \(\pi\) in C/C++ must be replaced by pi for Expression objects.) In the present case, the output of u_code and f_code is

x[0] + 2*x[1] + 1
-10*x[0] - 20*x[1] - 10

After having defined the mesh, the function space, and the boundary, we define the boundary value u_D as

u_D = Expression(u_code, degree=1)

Similarly, we define the right-hand side function as

f = Expression(f_code, degree=1)

FEniCS implementation¶

A solver for the nonlinear Poisson equation is as easy to implement as a solver for the linear Poisson equation. All we need to do is to state the formula for \(F\) and call solve(F == 0, u, bc) instead of solve(a == L, u, bc) as we did in the linear case. Here is a minimalistic code:

from fenics import *

def q(u):
    return 1 + u**2

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
u_D = Expression(u_code, degree=1)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

u = Function(V)
v = TestFunction(V)
f = Expression(f_code, degree=1)
F = q(u)*dot(grad(u), grad(v))*dx - f*v*dx

solve(F == 0, u, bc)

A complete version of this example program can be found in the file ft05_poisson_nonlinear.py.

The major difference from a linear problem is that the unknown function u in the variational form in the nonlinear case must be defined as a Function, not as a TrialFunction. In some sense this is a simplification from the linear case where we must define u first as a TrialFunction and then as a Function.

The solve function takes the nonlinear equations, derives symbolically the Jacobian matrix, and runs a Newton method to compute the solution.

# Warning: from fenics import * will import both `sym` and
# `q` from FEniCS. We therefore import FEniCS first and then
# overwrite these objects.
from fenics import *

def q(u):
    "Return nonlinear coefficient"
    return 1 + u**2

# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(q(u)*sym.diff(u, x), x) - sym.diff(q(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression(u_code, degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V)  # Note: not TrialFunction!
v = TestFunction(V)
f = Expression(f_code, degree=2)
F = q(u)*dot(grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc)

# Plot solution
plot(u)

# Compute maximum error at vertices. This computation illustrates
# an alternative to using compute_vertex_values as in poisson.py.
u_e = interpolate(u_D, V)
import numpy as np
error_max = np.abs(u_e.vector().array() - u.vector().array()).max()
print('error_max = ', error_max)

# Hold plot
interactive()

When we run the code, FEniCS reports on the progress of the Newton iterations. With \(2\cdot(8\times 8)\) cells, we reach convergence in eight iterations with a tolerance of \(10^{-9}\), and the error in the numerical solution is about \(10^{-16}\). These results bring evidence for a correct implementation. Thinking in terms of finite differences on a uniform mesh, \(\mathsf{P}_1\) elements mimic standard second-order differences, which compute the derivative of a linear or quadratic function exactly. Here, \(\nabla u\) is a constant vector, but then multiplied by \((1+u^2)\), which is a second-order polynomial in \(x\) and \(y\), which the divergence “difference operator” should compute exactly. We can therefore, even with \(\mathcal{P}_1\) elements, expect the manufactured \(u\) to be reproduced by the numerical method. With a nonlinearity like \(1+u^4\), this will not be the case, and we would need to verify convergence rates instead.

The current example shows how easy it is to solve a nonlinear problem in FEniCS. However, experts on the numerical solution of nonlinear PDEs know very well that automated procedures may fail for nonlinear problems, and that it is often necessary to have much better manual control of the solution process than what we have in the current case. We return to this problem in [Ref28] and show how we can implement tailored solution algorithms for nonlinear equations and also how we can steer the parameters in the automated Newton method used above.

Test problem¶

As a test example, we will model a clamped beam deformed under its own weight in 3D. This can be modeled by setting the right-hand side body force per unit volume to \(f=(0,0,-\varrho g)\) with \(\varrho\) the density of the beam and \(g\) the acceleration of gravity. The beam is box-shaped with length \(L\) and has a square cross section of width \(W\). We set \(u=u_{_\mathrm{D}} = (0,0,0)\) at the clamped end, \(x=0\). The rest of the boundary is traction free; that is, we set \(T = 0\).

FEniCS implementation¶

We first list the code and then comment upon the new constructions compared to the previous examples we have seen.

from fenics import *

# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u, title='Displacement', mode='displacement')

# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())

This example program can be found in the file ft06_elasticity.py.

Vector function spaces¶

The primary unknown is now a vector field \(u\) and not a scalar field, so we need to work with a vector function space:

V = VectorFunctionSpace(mesh, 'P', 1)

With u = Function(V) we get u as a vector-valued finite element function with three components for this 3D problem.

Constant vectors¶

For the boundary condition \(u=(0, 0, 0)\), we must set a vector value to zero, not just a scalar. Such a vector constant is specified as Constant((0, 0, 0)) in FEniCS. The corresponding 2D code would use Constant((0, 0)). Later in the code, we also need f as a vector and specify it as Constant((0, 0, rho*g)).

nabla_grad¶

The gradient and divergence operators now have a prefix nabla_. This is strictly not necessary in the present problem, but recommended in general for vector PDEs arising from continuum mechanics, if you interpret \(\nabla\) as a vector in the PDE notation; see the box about nabla_grad in the section Variational formulation.

Stress computation¶

As soon as the displacement u is computed, we can compute various stress measures. We will compute the von Mises stress defined as \(\sigma_M = \sqrt{\frac{3}{2}s:s}\) where \(s\) is the deviatoric stress tensor

There is a one-to-one mapping between these formulas and the FEniCS code:

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))

The von_Mises variable is now an expression that must be projected to a finite element space before we can visualize it:

V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')

Scaling¶

It is often advantageous to scale a problem as it reduces the need for setting physical parameters, and one obtains dimensionsless numbers that reflect the competition of parameters and physical effects. We develop the code for the original model with dimensions, and run the scaled problem by tweaking parameters appropriately. Scaling reduces the number of active parameters from 6 to 2 for the present application.

In Navier’s equation for \(u\), arising from inserting (37) and (38) into (36),

we insert coordinates made dimensionless by \(L\), and \(\bar u=u/U\), which results in the dimensionless governing equation

where \(\beta = 1 + \lambda/\mu\) is a dimensionless elasticity parameter and where

is a dimensionless variable reflecting the ratio of the load \(\varrho g\) and the shear stress term \(\mu\nabla^2 u\sim \mu U/L^2\) in the PDE.

One option for the scaling is to chose \(U\) such that \(\gamma\) is of unit size (\(U = \varrho gL^2/\mu\)). However, in elasticity, this leads to displacements of the size of the geometry, which makes plots look very strange. We therefore want the characteristic displacement to be a small fraction of the characteristic length of the geometry. This can be achieved by choosing \(U\) equal to the maximum deflection of a clamped beam, for which there actually exists a formula: \(U = \frac{3}{2}\varrho gL^2\delta^2/E\), where \(\delta = L/W\) is a parameter reflecting how slender the beam is, and \(E\) is the modulus of elasticity. Thus, the dimensionless parameter \(\delta\) is very important in the problem (as expected, since \(\delta\gg 1\) is what gives beam theory!). Taking \(E\) to be of the same order as \(\mu\), which is the case for many materials, we realize that \(\gamma \sim \delta^{-2}\) is an appropriate choice. Experimenting with the code to find a displacement that “looks right” in plots of the deformed geometry, points to \(\gamma = 0.4\delta^{-2}\) as our final choice of \(\gamma\).

The simulation code implements the problem with dimensions and physical parameters \(\lambda\), \(\mu\), \(\varrho\), \(g\), \(L\), and \(W\). However, we can easily reuse this code for a scaled problem: just set \(\mu = \varrho = L = 1\), \(W\) as \(W/L\) (\(\delta^{-1}\)), \(g=\gamma\), and \(\lambda=\beta\).

Plot of gravity-induced deflection in a clamped beam for the elasticity problem

Test problem 1: Channel flow¶

As a first test problem, we compute the flow between two infinite plates, so-called channel or Poiseuille flow. As we shall see, this problem has a known analytical solution. Let \(H\) be the distance between the plates and \(L\) the length of the channel. There are no body forces.

We may scale the problem first to get rid of seemingly independent physical parameters. The physics of this problem is governed by viscous effects only, in the direction perpendicular to the flow, so a time scale should be based on diffusion accross the channel: \(t_c = H^2/\nu\). We let \(U\), some characteristic inflow velocity, be the velocity scale and \(H\) the spatial scale. The pressure scale is taken as the characteristic shear stress, \(\mu U/H\), since this is a primary example of shear flow. Inserting \(\bar x = x/H\), \(\bar y = y/H\), \(\bar z = z/H\), \(\bar u =u/U\), \(\bar p = Hp/(\mu U)\), and \(\bar t = H^2/\nu\) in the equations results in the scaled Navier - Stokes equations (dropping bars after the scaling):

Here, \(\mathrm{Re} = \varrho UH/\mu\) is the Reynolds number. Because of the time and pressure scales, which are different from convection-dominated fluid flow, the Reynolds number is associated with the convective term and not the viscosity term.

The exact solution is derived by assuming \(u=(u_x(x,y,z),0,0)\), with the \(x\) axis pointing along the channel. Since \(\nabla\cdot u=0\), \(u\) cannot depend on \(x\). The physics of channel flow is also two-dimensional so we can omit the \(z\) coordinate (more precisely: \(\partial/\partial z=0\)). Inserting \(u=(u_x,0,0)\) in the (scaled) governing equations gives \(u_x''(y) = \partial p/\partial x\). Differentiating this equation with respect to \(x\) shows that \(\partial^2 p/\partial^2 x =0\) so \(\partial p/\partial x\) is a constant, here called \(-\beta\). This is the driving force of the flow and can be specified as a known parameter in the problem. Integrating \(u_x''(y)=-\beta\) over the width of the channel, \([0,1]\), and requiring \(u=(0, 0, 0)\) at the channel walls, results in \(u_x=\frac{1}{2}\beta y(1-y)\). The characteristic inlet velocity \(U\) can be taken as the maximum inflow at \(y=1/2\), implying \(\beta = 8\). The length of the channel, \(L/H\) in the scaled model, has no impact on the result, so for simplicity we just compute on the unit square. Mathematically, the pressure must be prescribed at a point, but since \(p\) does not depend on \(y\), we can set \(p\) to a known value, e.g. zero, along the outlet boundary \(x=1\). The result is \(p(x)=8(1-x)\) and \(u_x=4y(1-y)\).

The boundary conditions can be set as \(p=8\) at \(x=0\), \(p=0\) at \(x=1\) and \(u=(0, 0, 0)\) on the walls \(y=0,1\). This defines the pressure drop and should result in unit maximum velocity at the inlet and outlet and a parabolic velocity profile without further specifications. Note that it is only meaningful to solve the Navier - Stokes equations in 2D or 3D geometries, although the underlying mathematical problem collapses to two 1D problems, one for \(u_x(y)\) and one for \(p(x)\).

The scaled model is not so easy to simulate using a standard Navier - Stokes solver with dimensions. However, one can argue that the convection term is zero, so the Re coefficient in front of this term in the scaled PDEs is not important and can be set to unity. In that case, setting \(\varrho = \mu = 1\) in the original Navier - Stokes equations resembles the scaled model.

For a specific engineering problem one wants to simulate a specific fluid and set corresponding parameters. A general solver is therefore most naturally implemented with dimensions and the original physical parameters. However, scaling may greatly simplify numerical simulations. First of all, it shows that all fluids behave in the same way: it does not matter whether we have oil, gas, or water flowing between two plates, and it does not matter how fast the flow is (up to some criticial value of the Reynolds number where the flow becomes unstable and transitions to a complicated turbulent flow of totally different nature). This means that one simulation is enough to cover all types of channel flow! In other applications, scaling shows that it might be necessary to set just the fraction of some parameters (dimensionless numbers) rather than the parameters themselves. This simplifies exploring the input parameter space which is often the purpose of simulation. Frequently, the scaled problem is run by setting some of the input parameters with dimension to fixed values (often unity).

FEniCS implementation¶

Our previous examples have all started out with the creation of a mesh and then the definition of a FunctionSpace on the mesh. For the Navier - Stokes splitting scheme we will need to define two function spaces, one for the velocity and one for the pressure:

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

The first space V is a vector-valued function space for the velocity and the second space Q is a scalar-valued function space for the pressure. We use piecewise quadratic elements for the velocity and piecewise linear elements for the pressure. When creating a VectorFunctionSpace in FEniCS, the value-dimension (the length of the vectors) will be set equal to the geometric dimension of the finite element mesh. One can easily create vector-valued function spaces with other dimensions in FEniCS by adding the keyword parameter dim:

V = VectorFunctionSpace(mesh, 'P', 2, dim=10)

Since we have two different function spaces, we need to create two sets of trial and test functions:

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

As we have seen in previous examples, boundaries may be defined in FEniCS by defining Python functions that return True or False depending on whether a point should be considered part of the boundary, for example

def boundary(x, on_boundary):
    return near(x[0], 0)

This function defines the boundary to be all points with \(x\)-coordinate equal to (near) zero. The near function comes from FEniCS and performs a test with tolerance: abs(x[0] - 0) < 3E-16 so we do not run into rounding troubles. Alternatively, we may give the boundary definition as a string of C++ code, much like we have previously defined expressions such as u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2). The above definition of the boundary in terms of a Python function may thus be replaced by a simple C++ string:

boundary = 'near(x[0], 0)'

This has the advantage of moving the computation of which nodes belong to the boundary from Python to C++, which improves the efficiency of the program.

For the current example, we will set three different boundary conditions. First, we will set \(u = 0\) at the walls of the channel; that is, at \(y = 0\) and \(y = 1\). Second, we will set \(p = 8\) at the inflow (\(x = 0\)) and, finally, \(p = 0\) at the outflow (\(x = 1\)). This will result in a pressure gradient that will accelerate the flow from the initial state with zero velocity. These boundary conditions may be defined as follows:

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 1)'
walls    = 'near(x[1], 0) || near(x[1], 1)'

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

At the end, we collect the boundary conditions for the velocity and pressure in Python lists so we can easily access them in the following computation.

We now move on to the definition of the variational forms. There are three variational problems to be defined, one for each step in the IPCS scheme. Let us look at the definition of the first variational problem. We start with some constants:

U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

The next step is to set up the variational form for the first step (48) in the solution process. Since the variational problem contains a mix of known and unknown quantities we will use the following naming convention: u is the unknown (mathematically \(u^{n+1}\)) as a trial function in the variational form, u_ is the most recently computed approximation (\(u^{n+1}\) available as a Function object), u_n is \(u^n\), and the same convention goes for p, p_ (\(p^{n+1}\)), and p_n (\(p^n\)).

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
     rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

Note that we take advantage of the Python programming language to define our own operators sigma and epsilon. Using Python this way makes it easy to extend the mathematical language of FEniCS with special operators and constitutive laws.

Also note that FEniCS can sort out the bilinear form \(a(u,v)\) and linear form \(L(v)\) forms by the lhs and rhs functions. This is particularly convenient in longer and more complicated variational forms.

The splitting scheme requires the solution of a sequence of three variational problems in each time step. We have previously used the built-in FEniCS function solve to solve variational problems. Under the hood, when a user calls solve(a == L, u, bc), FEniCS will perform the following steps:

A = assemble(A)
b = assemble(L)
bc.apply(A, b)
solve(A, u.vector(), b)

In the last step, FEniCS uses the overloaded solve function to solve the linear system AU = b where U is the vector of degrees of freedom for the function \(u(x) = \sum_{j=1} U_j \phi_j(x)\).

In our implementation of the splitting scheme, we will make use of these low-level commands to first assemble and then call solve. This has the advantage that we may control when we assemble and when we solve the linear system. In particular, since the matrices for the three variational problems are all time-independent, it makes sense to assemble them once and for all outside of the time-stepping loop:

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

Within the time-stepping loop, we may then assemble only the right-hand side vectors, apply boundary conditions, and call the solve function as here for the first of the three steps:

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1)

Notice the Python list comprehension [bc.apply(b1) for bc in bcu] which iterates over all bc in the list bcu. This is a convenient and compact way to construct a loop that applies all boundary conditions in a single line. Also, the code works if we add more Dirichlet boundary conditions in the future. Note that the boundary conditions only need to be applied to the right-hand side vectors as they have already been applied to the matrices (not shown).

Finally, let us look at an important detail in how we use parameters such as the time step dt in the definition of our variational problems. Since we might want to change these later, for example if we want to experiment with smaller or larger time steps, we wrap these using a FEniCS Constant:

k = Constant(dt)

The assembly of matrices and vectors in FEniCS is based on code generation. This means that whenever we change a variational problem, FEniCS will have to generate new code, which may take a little time. New code will also be generated and compiled when a float value for the time step is changed. By wrapping this parameter using Constant, FEniCS will treat the parameter as a generic constant and not as a specific numerical value, which prevents repeated code generation. In the case of the time step, we choose a new name k instead of dt for the Constant since we also want to use the variable dt as a Python float as part of the time-stepping.

The complete code for simulating 2D channel flow with FEniCS can be found in the file ft07_navier_stokes_channel.py.

Verification¶

We compute the error at the nodes as we have done before to verify that our implementation is correct. Our Navier - Stokes solver computes the solution to the time-dependent incompressible Navier - Stokes equations, starting from the initial condition \(u = (0, 0)\). We have not specified the initial condition explicitly in our solver which means that FEniCS will initialize all variables, in particular the previous and current velocities u_n and u_, to zero. Since the exact solution is quadratic, we expect the solution to be exact to within machine precision at the nodes at infinite time. For our implementation, the error quickly approaches zero and is approximately \(10^{-6}\) at time \(T = 10\).

Plot of the velocity profile at the final time for the Navier - Stokes channel flow example

Test problem 2: Flow past a cylinder¶

We now turn our attention to a more challenging problem: flow past a circular cylinder. The geometry and parameters are taken from problem DFG 2D-2 in the FEATFLOW/1995-DFG benchmark suite and is illustrated in Figure Geometry for the flow past a cylinder test problem. Notice the slightly perturbed and unsymmetric geometry. The kinematic viscosity is given by \(\nu = 0.001 = \mu/\varrho\) and the inflow velocity profile is specified as

which has a maximum magnitude of \(1.5\) at \(y = 0.41/2\). We do not use any scaling for this problem since all exact parameters are known.

Geometry for the flow past a cylinder test problem. Notice the slightly perturbed and unsymmetric geometry

FEniCS implementation¶

So far all our domains have been simple shapes such as a unit square or a rectangular box. A number of such simple meshes may be created using the built-in mesh classes in FEniCS (UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh, IntervalMesh, RectangleMesh, BoxMesh). FEniCS supports the creation of more complex meshes via a technique called constructive solid geometry (CSG), which lets us define geometries in terms of simple shapes (primitives) and set operations: union, intersection, and set difference. The set operations are encoded in FEniCS using the operators + (union), * (intersection), and - (set difference). To access the CSG functionality in FEniCS, one must import the FEniCS module mshr which provides the extended meshing functionality of FEniCS.

The geometry for the cylinder flow test problem can be defined easily by first defining the rectangular channel and then subtracting the circle:

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder

We may then create the mesh by calling the function generate_mesh:

mesh = generate_mesh(domain, 64)

Here the argument 64 indicates that we want to resolve the geometry with 64 cells across its diameter (the channel length).

To solve the cylinder test problem, we only need to make a few minor changes to the code we wrote for the channel flow test case. Besides defining the new mesh, the only change we need to make is to modify the boundary conditions and the time step size. The boundaries are specified as follows:

inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

The last line may seem cryptic before you catch the idea: we want to pick out all boundary points (on_boundary) that also lie within the 2D domain \([0.1,0.3]\times [0.1,0.3]\), see Figure Geometry for the flow past a cylinder test problem. Notice the slightly perturbed and unsymmetric geometry. The only possible points are then the points on the circular boundary!

In addition to these essential changes, we will make a number of small changes to improve our solver. First, since we need to choose a relatively small time step to compute the solution (a time step that is too large will make the solution blow up) we add a progress bar so that we can follow the progress of our computation. This can be done as follows:

progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0.0
for n in range(num_steps):

    # Update current time
    t += dt

    # Place computation here

    # Update progress bar
    progress.update(t / T)

Since the system(s) of linear equations are significantly larger than for the simple channel flow test problem, we choose to use an iterative method instead of the default direct (sparse) solver used by FEniCS when calling solve. Efficient solution of linear systems arising from the discretization of PDEs requires the choice of both a good iterative (Krylov subspace) method and a good preconditioner. For this problem, we will simply use the biconjugate gradient stabilized method (BiCGSTAB) and the conjugate gradient method. This can be done by adding the keywords bicgstab or cg in the call to solve. We also specify suitable preconditioners to speed up the computations:

solve(A1, u1.vector(), b1, 'bicgstab', 'hypre_amg')
solve(A2, p1.vector(), b2, 'bicgstab', 'hypre_amg')
solve(A3, u1.vector(), b3, 'cg', 'sor')

Finally, to be able to postprocess the computed solution in ParaView, we store the solution to a file in each time step. We have previously created files with the suffix .pvd for this purpose. In the example program ft04_heat_gaussian.py, we first created a file named heat_gaussian/solution.pvd and then saved the solution in each time step using

vtkfile << (u, t)

For the present example, we will instead choose to save the solution to XDMF format. This file format works similarly to the .pvd files we have seen earlier but has several advantages. First, the storage is much more efficient, both in terms of speed and file sizes. Second, .xdmf files work in parallell, both for writing and reading (postprocessing). Much like .pvd files, the actual data will not be stored in the .xdmf file itself, but will instead be stored in a (single) separate data file with the suffix .hdf5 which is an advanced file format designed for high-performance computing. We create the XDMF files as follows:

xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

In each time step, we may then store the velocity and pressure by

xdmffile_u.write(u, t)
xdmffile_p.write(p, t)

We also store the solution using a FEniCS TimeSeries. This allows us to store the solution not for visualization, but for later reuse in a computation as we will see in the next section. Using a TimeSeries it is easy and efficient to read in solutions from certain points in time during a simulation. The TimeSeries class also uses the HDF5 file format for efficient storage and access to data.

Figures Plot of the velocity for the cylinder test problem at final time and Plot of the pressure for the cylinder test problem at final time show the velocity and pressure at final time visualized in ParaView. For the visualization of the velocity, we have used the Glyph filter to visualize the vector velocity field. For the visualization of the pressure, we have used the Warp By Scalar filter.

Plot of the velocity for the cylinder test problem at final time

Plot of the pressure for the cylinder test problem at final time

The complete code for the cylinder test problem looks as follows:

from fenics import *
from mshr import *
import numpy as np

T = 5.0            # final time
num_steps = 5000   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    plot(u_, title='Velocity')
    plot(p_, title='Pressure')

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    progress.update(t / T)
    print('u max:', u_.vector().array().max())

# Hold plot
interactive()

This program can be found in the file ft08_navier_stokes_cylinder.py. The reader should be advised that this example program is considerably more demanding than our previous examples in terms of CPU time and memory, but it should be possible to run the program on a reasonably modern laptop.

Setting initial conditions for mixed systems¶

In our example, we did not need to worry about setting an initial condition, since we start with \(u_1 = u_2 = u_3 = 0\). This happens automatically in the code when we set u_n = Function(V). This creates a Function for the whole system and all degrees of freedom are set to zero.

If we want to set initial conditions for the components of the system separately, the easiest solution is to define the initial conditions as a vector-valued Expression and then project (or interpolate) this to the Function representing the whole system. For example,

u_0 = Expression(('sin(x[0])', 'cos(x[0]*x[1])', 'exp(x[1])'), degree=1)
u_n = project(u_0, V)

This defines \(u_1\), \(u_2\), and \(u_2\) to be the projections of \(\sin x\), \(\cos (xy)\), and \(\exp(y)\), respectively.

Setting boundary conditions for mixed systems¶

In our example, we also did not need to worry about setting boundary conditions since we used a natural Neumann condition. If we want to set Dirichlet conditions for individual components of the system, this can be done as usual by the class DirichletBC, but we must specify for which subsystem we set the boundary condition. For example, to specify that \(u_2\) should be equal to \(xy\) on the boundary defined by boundary, we do

u_D = Expression('x[0]*x[1]', degree=1)
bc = DirichletBC(V.sub(1), u_D, boundary)

The object bc or a list of such objects containing different boundary conditions, can then be passed to the solve function as usual. Note that numbering starts at \(0\) in FEniCS so the subspace corresponding to \(u_2\) is V.sub(1).

Accessing components of mixed systems¶

If u is a Function defined on a mixed function space in FEniCS, there are several ways in which u can be split into components. Above we already saw an example of the first of these:

u_1, u_2, u_3 = split(u)

This extracts the components of u as symbols that can be used in a variational problem. The above statement is in fact equivalent to

u_1 = u[0]
u_2 = u[1]
u_3 = u[2]

Note that u[0] is not really a Function object, but merely a symbolic expression, just like grad(u) in FEniCS is a symbolic expression and not a Function representing the gradient. This means that u_1, u_2, u_3 can be used in a variational problem, but cannot be used for plotting or postprocessing.

To access the components of u for plotting and saving the solution to file, we need to use a different variant of the split function:

u_1_, u_2_, u_3_ = u.split()

This returns three subfunctions as actual objects with access to the common underlying data stored in u, which makes plotting and saving to file possible. Alternatively, we can do

u_1_, u_2_, u_3_ = u.split(deepcopy=True)

which will create u_1_, u_2_, and u_3_ as stand-alone Function objects, each holding a copy of the subfunction data extracted from u. This is useful in many situations but is not necessary for plotting and saving solutions to file.