Stokes equation

Find $u \in [H^1_D]^2$ and $p \in L_2$ such that

$$ \DeclareMathOperator{\Div}{div} \begin{array}{ccccll} \int \nabla u : \nabla v & + & \int \Div v \, p & = & \int f v & \forall \, v \\ \int \Div u \, q & & & = & 0 & \forall \, q \end{array} $$

Define channel geometry and mesh it:

In [1]:
from ngsolve import *
import netgen.gui
%gui tk

from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl")
mesh = Mesh( geo.GenerateMesh(maxh=0.05))
mesh.Curve(3)
Draw (mesh)

bcs =  ('wall', 'outlet', 'wall', 'inlet')

Use Taylor Hood finite element pairing: Continuous $P^2$ elements for velocity, and continuous $P^1$ for pressure:

In [2]:
V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=1)
X = FESpace([V,V,Q])

Setup bilinear-form for Stokes. We give names for all scalar field components. The divergence is constructed from partial derivatives of the velocity components.

In [3]:
ux,uy,p = X.TrialFunction()
vx,vy,q = X.TestFunction()

div_u = grad(ux)[0]+grad(uy)[1]
div_v = grad(vx)[0]+grad(vy)[1]

a = BilinearForm(X)
a += SymbolicBFI(grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p)
a.Assemble()

Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

In [4]:
gfu = GridFunction(X)
uin = 1.5*4*y*(0.41-y)/(0.41*0.41)
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
velocity = CoefficientFunction(gfu.components[0:2])
Draw(velocity, mesh, "vel")
SetVisualization(max=2)

Solve equation:

In [5]:
res = gfu.vec.CreateVector()
res.data = -a.mat * gfu.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
gfu.vec.data += inv * res
Redraw()

Now we define a Stokes setup function to test different spaces:

In [6]:
def SolveStokes(X):
    ux,uy,p = X.TrialFunction()
    vx,vy,q = X.TestFunction()
    div_u = grad(ux)[0]+grad(uy)[1]
    div_v = grad(vx)[0]+grad(vy)[1]
    a = BilinearForm(X)
    a += SymbolicBFI(grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p)
    a.Assemble()
    gfu = GridFunction(X)
    uin = 1.5*4*y*(0.41-y)/(0.41*0.41)
    gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
    res = gfu.vec.CreateVector()
    res.data = -a.mat * gfu.vec
    inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
    gfu.vec.data += inv * res
    return gfu

Higher order Taylor-Hood elements:

In [7]:
V = H1(mesh, order=4, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=3)
X = FESpace([V,V,Q])

gfu = SolveStokes(X)

velocity = CoefficientFunction(gfu.components[0:2])
Draw(velocity, mesh, "vel")
SetVisualization(max=2)

Discontinuous pressure elements:

In [8]:
V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = FESpace([V,V,Q])

try:
    gfu = SolveStokes(X)
except RuntimeError as a:
    print(a)

velocity = CoefficientFunction(gfu.components[0:2])
Draw(velocity, mesh, "vel")

V.ndof = 1694 , Q.ndof = 2370
UmfpackInverse: Numeric factorization failed.

$P^{2,+} \times P^{1,dc}$ elements:

In [9]:
V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
V.order[TRIG]=3
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
Q = L2(mesh, order=1)
X = FESpace([V,V,Q])

gfu = SolveStokes(X)

velocity = CoefficientFunction(gfu.components[0:2])
Draw(velocity, mesh, "vel")

V.ndof = 2484 , Q.ndof = 2370

the mini element:

In [10]:
V = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
V.order[TRIG]=3
Q = H1(mesh, order=1)
X = FESpace([V,V,Q])

gfu = SolveStokes(X)

velocity = CoefficientFunction(gfu.components[0:2])
Draw(velocity, mesh, "vel")

A newer feature: A vector-valued $H^1$-space: Less to type and more possibilities to explore structure and optimize.

In [11]:
V = FESpace("VectorH1", mesh, order=3, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=2)
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

a = BilinearForm(X)
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p)
a.Assemble()

gfu = GridFunction(X)
uin = 1.5*4*y*(0.41-y)/(0.41*0.41)
gfu.components[0].components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

res = gfu.vec.CreateVector()
res.data = -a.mat * gfu.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
gfu.vec.data += inv * res
Draw(gfu.components[0], mesh, "vel")
In [12]: