21. Stokes equations with stabilized first order elements

This demo illustrates how to:

  • Solve the Stokes problem using a stabilized method

The velocity \(u\) and \(p\) will look as follows:

../../../../_images/plot_u3.png ../../../../_images/plot_p1.png

21.1. Equation and problem definition

Strong formulation of the problem:

\[\begin{split}- \nabla \cdot (\nabla u + p I) &= f \quad {\rm in} \ \Omega, \\ \nabla \cdot u &= 0 \quad {\rm in} \ \Omega. \\\end{split}\]

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:

\[\begin{split}u &= u_0 \quad {\rm on} \ \Gamma_{D}, \\ \nabla u \cdot n + p n &= g \quad {\rm on} \ \Gamma_{N}. \\\end{split}\]

The Stokes equations can easily be 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

\[a((u, p), (v, q)) = L((v, q))\]

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

\[\begin{split}a((u, p), (v, q)) &= \int_{\Omega} \nabla u \cdot \nabla v - \nabla \cdot v \ p + \nabla \cdot u \ q \, {\rm d} x, \\ L((v, q)) &= \int_{\Omega} f \cdot v \, {\rm d} x \, {\rm d} s. \\\end{split}\]

The space \(W\) should be a mixed (product) function space: \(W = V \times Q\) such that \(u \in V\) and \(q \in Q\). In this demo we will use first order elements in both velocity and pressure, which will lead to stability problems. Therefore we use a stabilized variational formulation. It reads as follows: find \((u, p) \in W\) such that

\[a((u, p), (v, q)) = L((v, q))\]

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

\[\begin{split}a((u, p), (v, q)) &= \int_{\Omega} \nabla u \cdot \nabla v - \nabla \cdot v \ p + \nabla \cdot u \ q + \delta \nabla q \cdot \nabla p \, {\rm d} x, \\ L((v, q)) &= \int_{\Omega} f \cdot v \, {\rm d} x + \delta \nabla q \cdot f \, {\rm d} s.\end{split}\]

where \(\delta = \beta h^2\) and \(\beta\) is some number and \(h\) is the mesh cell size.

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

  • \((\Omega = [0,1]\times[0, 1]) \backslash {\rm dolphin}\)
  • \(\Gamma_D =\)
  • \(\Gamma_N =\)
  • \(u_0 = (- \sin(\pi x_1), 0.0) \, {\rm for} \, x_0 = 1 \, {\rm and} \, u_0 = (0.0, 0.0)\) otherwise
  • \(f = (0.0, 0.0)\)
  • \(g = (0.0, 0.0)\)

21.2. Implementation

First, the dolfin module is imported:

from dolfin import *

In this example, different boundary conditions are prescribed on different parts of the boundaries. This information must be made available to the solver. One way of doing this, is to tag the different sub-regions with different (integer) labels. DOLFIN provides a class MeshFunction which is useful for these types of operations: instances of this class represent functions over mesh entities (such as over cells or over facets). Mesh and mesh functions can be read from file in the following way:

# Load mesh and subdomains
mesh = Mesh("dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "dolfin_fine_subdomains.xml.gz")

Next, we define a MixedFunctionSpace composed of a VectorFunctionSpace and a FunctionSpace, both of continuous piecewise linears.

# Define function spaces
scalar = FunctionSpace(mesh, "CG", 1)
vector = VectorFunctionSpace(mesh, "CG", 1)
system = vector * scalar

Now that we have our mixed function space and marked subdomains defining the boundaries, we create functions for the boundary conditions and define boundary conditions:

# Create functions for boundary conditions
noslip = Constant((0, 0))
inflow = Expression(("-sin(x[1]*pi)", "0"))
zero   = Constant(0)

# No-slip boundary condition for velocity
bc0 = DirichletBC(system.sub(0), noslip, sub_domains, 0)

# Inflow boundary condition for velocity
bc1 = DirichletBC(system.sub(0), inflow, sub_domains, 1)

# Boundary condition for pressure at outflow
bc2 = DirichletBC(system.sub(1), zero, sub_domains, 2)

# Collect boundary conditions
bcs = [bc0, bc1, bc2]

Here, we have given four arguments to DirichletBC. The first specifies the FunctionSpace. Since we have a MixedFunctionSpace, we write system.sub(0) for the VectorFunctionSpace, and system.sub(1) for the FunctionSpace. The second argument specifies the value on the Dirichlet boundary. The two last arguments specify the marking of the subdomains; sub_domains contains the subdomain markers and the number given as the last argument is the subdomain index.

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

# Define variational problem
(v, q) = TestFunctions(system)
(u, p) = TrialFunctions(system)
f = Constant((0, 0))
h = CellSize(mesh)
beta  = 0.2
delta = beta*h*h
a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u) + \
    delta*inner(grad(q), grad(p)))*dx
L = inner(v + delta*grad(q), f)*dx

To compute the solution we use the bilinear and linear forms, and the boundary condition, but we also need to create a Function to store the solution(s). The (full) solution will be stored in w, which we initialize using the MixedFunctionSpace system. The actual computation is performed by calling solve with the arguments a, L and bcs. The separate components u and p of the solution can be extracted by calling the split function.

# Compute solution
w = Function(system)
solve(a == L, w, bcs)
u, p = w.split()

Finally, we can store to file and plot the solutions.

# 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()

21.3. Complete code


from dolfin import *

# Load mesh and subdomains
mesh = Mesh("../dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "../dolfin_fine_subdomains.xml.gz");

# Define function spaces
scalar = FunctionSpace(mesh, "CG", 1)
vector = VectorFunctionSpace(mesh, "CG", 1)
system = vector * scalar

# Create functions for boundary conditions
noslip = Constant((0, 0))
inflow = Expression(("-sin(x[1]*pi)", "0"))
zero   = Constant(0)

# No-slip boundary condition for velocity
bc0 = DirichletBC(system.sub(0), noslip, sub_domains, 0)

# Inflow boundary condition for velocity
bc1 = DirichletBC(system.sub(0), inflow, sub_domains, 1)

# Boundary condition for pressure at outflow
bc2 = DirichletBC(system.sub(1), zero, sub_domains, 2)

# Collect boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
(v, q) = TestFunctions(system)
(u, p) = TrialFunctions(system)
f = Constant((0, 0))
h = CellSize(mesh)
beta  = 0.2
delta = beta*h*h
a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u) + \
    delta*inner(grad(q), grad(p)))*dx
L = inner(v + delta*grad(q), f)*dx

# Compute solution
w = Function(system)
solve(a == L, w, bcs)
u, p = w.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()