# 18. Stokes equations¶

This demo is implemented in a single Python file,
`demo_stokes-iterative.py`

, which contains both the
variational forms and the solver.

This demo illustrates how to:

- Maintain symmetry when assembling a system of symmetric equations with essential (Dirichlet) boundary conditions
- Use a iterative solver explicitly for solving a linear system of equations
- Define a preconditioner explicitly using a form

## 18.1. Strong form of the Stokes equations¶

The incompressible Stokes equations in strong form read: for a domain \(\Omega \subset \mathbb{R}^n\), find the velocity \(u\) and the pressure \(p\) satisfying

Note

The sign of the pressure has been flipped from the classical definition. This is done in order to have a symmetric (but not positive-definite) system of equations rather than a non-symmetric (but positive-definite) system of equations.

A typical set of boundary conditions on the boundary \(\partial \Omega = \Gamma_{D} \cup \Gamma_{N}\) can be:

## 18.2. Weak form of the Stokes equations¶

The Stokes equations can easily formulated in a mixed variational form; that is, a form where the two variables, the velocity and the pressure, are approximated simultaneously. Using the abstract framework, we have the problem: find \((u, p) \in W\) such that

for all \((v, q) \in W\) where

The space W should be a mixed (product) function space: \(W = V \times Q\) such that \(u \in V\) and \(q \in Q\).

## 18.3. Preconditioning of the linear system of equations¶

For the resulting linear system of equations, the following form defines a suitable preconditioner:

## 18.4. Domain and boundary conditions¶

In this demo, we shall consider the following definitions of the input functions, the domain, and the boundaries:

- \(\Omega = [0,1]^3\) (a unit cube)
- \(\Omega_D = \{(x_0, x_1, x_2) \, | \, x_0 = 0 \, \text{or} \, x_0 = 1 \, \text{or} \, x_1 = 0 \, \text{or} \, x_1 = 1 \}\)
- \(u_0 = (- \sin(\pi x_1), 0.0, 0.0)\) for \(x_0 = 1\) and \(u_0 = (0.0, 0.0, 0.0)\) otherwise
- \(f = (0.0, 0.0, 0.0)\)
- \(g = (0.0, 0.0, 0.0)\)

## 18.5. Implementation¶

This description goes through the implementation (in
`demo_stokes-iterative.py`

) of a solver for the above
described Stokes equations. Some of the standard steps will be
described in less detail, so before reading this, we suggest that you
are familiarize with the Poisson demo (for the very basics) and the
Mixed Poisson demo (for how to deal with
mixed function spaces). Also, the Navier–Stokes demo illustrates how to use
iterative solvers in a more implicit manner (typically only suitable
for positive-definite systems of equations).

The Stokes equations as formulated above result in a system of linear equations that is not positive-definite. Standard iterative linear solvers typically fail to converge for such systems. Some care must therefore be taken in preconditioning the systems of equations. Moreover, not all of the linear algebra backends support this. We therefore start by checking that either “PETSc” or “Epetra” (from Trilinos) is available:

```
from dolfin import *
# Test for PETSc or Epetra
if not has_linear_algebra_backend("PETSc") and not has_linear_algebra_backend("Epetra"):
info("DOLFIN has not been configured with Trilinos or PETSc. Exiting.")
exit()
if not has_krylov_solver_preconditioner("amg"):
info("Sorry, this demo is only available when DOLFIN is compiled with AMG "
"preconditioner, Hypre or ML.");
exit()
```

Next, we define the mesh (a `UnitCubeMesh`

) and a `MixedFunctionSpace`

composed of a
`VectorFunctionSpace`

of continuous
piecewise quadratics and a `FunctionSpace`

of continuous
piecewise linears. (This mixed finite element space is known as the
Taylor–Hood elements and is a stable, standard element pair for the
Stokes equations.)

```
# Load mesh
mesh = UnitCubeMesh(16, 16, 16)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q
```

Next, we define the boundary conditions.

```
# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top_bottom(x, on_boundary):
return x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS
# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0", "0.0"))
bc1 = DirichletBC(W.sub(0), inflow, right)
# Boundary condition for pressure at outflow
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, left)
# Collect boundary conditions
bcs = [bc0, bc1, bc2]
```

The bilinear and linear forms corresponding to the weak mixed formulation of the Stokes equations are defined as follows:

```
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
```

We can now use the same `TrialFunctions`

and
`TestFunctions`

to
define the preconditioner matrix. We first define the form
corresponding to the expression for the preconditioner (given in the
initial description above):

```
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
```

Next, we want to assemble the matrix corresponding to the bilinear
form and the vector corresponding to the linear form of the Stokes
equations. Moreover, we want to apply the specified boundary
conditions to the linear system. However, `assembling`

the matrix and vector and applying a
`DirichletBC`

separately
will possibly result in a non-symmetric system of equations. Instead,
we can use the `assemble_system`

function to assemble both the
matrix `A`

, the vector `bb`

, and apply the boundary conditions
`bcs`

in a symmetric fashion:

```
# Assemble system
A, bb = assemble_system(a, L, bcs)
```

We do the same for the preconditioner matrix `P`

using the linear
form `L`

as a dummy form:

```
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
```

Next, we specify the iterative solver we want to use, in this case a
`KrylovSolver`

. We associate the
left-hand side matrix `A`

and the preconditioner matrix `P`

with
the solver by calling `solver.set_operators`

.

```
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver("tfqmr", "amg")
# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
```

We are now almost ready to solve the linear system of equations. It
remains to specify a `Vector`

for
storing the result. For easy manipulation later, we can define a
`Function`

and use the
vector associated with this Function. The call to
`solver.solve`

then looks as
follows

```
# Solve
U = Function(W)
solver.solve(U.vector(), bb)
```

Finally, we can play with the result in different ways:

```
# Get sub-functions
u, p = U.split()
# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
# Plot solution
plot(u)
plot(p)
interactive()
```

## 18.6. Complete code¶

```
from dolfin import *
# Test for PETSc or Epetra
if not has_linear_algebra_backend("PETSc") and not has_linear_algebra_backend("Epetra"):
info("DOLFIN has not been configured with Trilinos or PETSc. Exiting.")
exit()
if not has_krylov_solver_preconditioner("amg"):
info("Sorry, this demo is only available when DOLFIN is compiled with AMG "
"preconditioner, Hypre or ML.");
exit()
# Load mesh
mesh = UnitCubeMesh(16, 16, 16)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q
# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top_bottom(x, on_boundary):
return x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS
# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0", "0.0"))
bc1 = DirichletBC(W.sub(0), inflow, right)
# Boundary condition for pressure at outflow
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, left)
# Collect boundary conditions
bcs = [bc0, bc1, bc2]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
# Assemble system
A, bb = assemble_system(a, L, bcs)
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver("tfqmr", "amg")
# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
# Solve
U = Function(W)
solver.solve(U.vector(), bb)
# Get sub-functions
u, p = U.split()
# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
# Plot solution
plot(u)
plot(p)
interactive()
```