This page was generated from unit-2.6-stokes/stokes.ipynb.

2.6 Stokes equation

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

\[\begin{split}\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}\end{split}\]

Define channel geometry and mesh it:

from ngsolve import *
from ngsolve.webgui import Draw

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)

BaseWebGuiScene

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

V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=1)
X = 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.

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 += (grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p) * dx
a.Assemble()

<ngsolve.comp.BilinearForm at 0x7f41781a4fb0>

Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

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:

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(velocity, mesh, "vel")

BaseWebGuiScene

Testing different velocity-pressure pairs

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

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 += (grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p)*dx
    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


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

    return gfu

Higher order Taylor-Hood elements:

V = H1(mesh, order=4, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=3)
X = V*V*Q

gfu = SolveStokes(X)

With discontinuous pressure elements P2-P1 is unstable:

V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = V*V*Q

gfu = SolveStokes(X)

V.ndof = 1656 , Q.ndof = 2322

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
<ipython-input-1-5ecc0da3557f> in <module>
      4 X = V*V*Q
      5
----> 6 gfu = SolveStokes(X)

<ipython-input-1-abe06cca85c3> in SolveStokes(X)
     12     res = gfu.vec.CreateVector()
     13     res.data = -a.mat * gfu.vec
---> 14     inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
     15     gfu.vec.data += inv * res
     16

NgException: UmfpackInverse: Numeric factorization failed.

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

V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
V.SetOrder(TRIG,3)
V.Update()
Q = L2(mesh, order=1)
X = V*V*Q
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)

gfu = SolveStokes(X)

V.ndof = 2430 , Q.ndof = 2322

the mini element:

V = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
V.SetOrder(TRIG,3)
V.Update()
Q = H1(mesh, order=1)
X = V*V*Q

gfu = SolveStokes(X)

VectorH1

A vector-valued \(H^1\)-space: Less to type and more possibilities to explore structure and optimize.

V = VectorH1(mesh, order=2, dirichlet="wall|inlet|cyl")
V.SetOrder(TRIG,3)
V.Update()
Q = L2(mesh, order=1)
X = V*Q

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

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

gfu = GridFunction(X)
uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
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
Draw(gfu.components[0], mesh, "vel")
Draw(Norm(gfu.components[0]), mesh, "|vel|")
SetVisualization(max=2)

Stokes as a block-system

We can now define separate bilinear-form and matrices for A and B, and combine them to a block-system:

V = VectorH1(mesh, order=3, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=2)

u,v = V.TnT()
p,q = Q.TnT()

a = BilinearForm(V)
a += InnerProduct(Grad(u),Grad(v))*dx

b = BilinearForm(trialspace=V, testspace=Q)
b += div(u)*q*dx

a.Assemble()
b.Assemble()

<ngsolve.comp.BilinearForm at 0x7f413eb2db70>

Needed as preconditioner for the pressure:

mp = BilinearForm(Q)
mp += p*q*dx
mp.Assemble()

<ngsolve.comp.BilinearForm at 0x7f413eaecf30>

Two right hand sides for the two spaces:

f = LinearForm(V)
f += CoefficientFunction((0,x-0.5)) * v * dx
f.Assemble()

g = LinearForm(Q)
g.Assemble()

<ngsolve.comp.LinearForm at 0x7f413f70cfb0>

Two GridFunctions for velocity and pressure:

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")
uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.Set(uin, definedon=mesh.Boundaries("inlet"))

Combine everything to a block-system. BlockMatrix and BlockVector store references to the original matrices and vectors, no new large matrices are allocated. The same for the transpose matrix b.mat.T. It stores a wrapper for the original matrix, and replaces the call of the Mult function by MultTrans.

K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, None] ] )
C = BlockMatrix( [ [a.mat.Inverse(V.FreeDofs()), None], [None, mp.mat.Inverse()] ] )

rhs = BlockVector ( [f.vec, g.vec] )
sol = BlockVector( [gfu.vec, gfp.vec] )

solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, initialize=False)

LinearSolver iteration 1, residual = 4.337069360146072
LinearSolver iteration 2, residual = 2.264713596298657
LinearSolver iteration 3, residual = 1.892551008538201
LinearSolver iteration 4, residual = 1.5278298950439118
LinearSolver iteration 5, residual = 1.4774123933006966
LinearSolver iteration 6, residual = 1.2748646703687811
LinearSolver iteration 7, residual = 1.227350244839074
LinearSolver iteration 8, residual = 1.0796748148348154
LinearSolver iteration 9, residual = 0.9927213467063855
LinearSolver iteration 10, residual = 0.8697728301343073
LinearSolver iteration 11, residual = 0.816271182836896
LinearSolver iteration 12, residual = 0.7259673235404897
LinearSolver iteration 13, residual = 0.7046418006753331
LinearSolver iteration 14, residual = 0.6436420998832291
LinearSolver iteration 15, residual = 0.6266519183176467
LinearSolver iteration 16, residual = 0.5916941058249604
LinearSolver iteration 17, residual = 0.5778872637584644
LinearSolver iteration 18, residual = 0.5491983683896285
LinearSolver iteration 19, residual = 0.5321339542232276
LinearSolver iteration 20, residual = 0.5051520176131404
LinearSolver iteration 21, residual = 0.49388035273938846
LinearSolver iteration 22, residual = 0.4774977661080504
LinearSolver iteration 23, residual = 0.46005121716109815
LinearSolver iteration 24, residual = 0.43657068469684934
LinearSolver iteration 25, residual = 0.41356565003754386
LinearSolver iteration 26, residual = 0.3795032471424406
LinearSolver iteration 27, residual = 0.3467761879111833
LinearSolver iteration 28, residual = 0.30529859256434966
LinearSolver iteration 29, residual = 0.25676736950680873
LinearSolver iteration 30, residual = 0.1935837308698832
LinearSolver iteration 31, residual = 0.14865328514139542
LinearSolver iteration 32, residual = 0.09766275665495895
LinearSolver iteration 33, residual = 0.0825927053478491
LinearSolver iteration 34, residual = 0.04796176411133548
LinearSolver iteration 35, residual = 0.046561374422869194
LinearSolver iteration 36, residual = 0.024448642105259134
LinearSolver iteration 37, residual = 0.024187420660011078
LinearSolver iteration 38, residual = 0.013787348456802321
LinearSolver iteration 39, residual = 0.013782015309911428
LinearSolver iteration 40, residual = 0.009231401649420294
LinearSolver iteration 41, residual = 0.009136165052020364
LinearSolver iteration 42, residual = 0.006992878478949685
LinearSolver iteration 43, residual = 0.006933444411962679
LinearSolver iteration 44, residual = 0.004323990920530351
LinearSolver iteration 45, residual = 0.004307641889957301
LinearSolver iteration 46, residual = 0.002585414701387047
LinearSolver iteration 47, residual = 0.002579809247615046
LinearSolver iteration 48, residual = 0.0016354435576020616
LinearSolver iteration 49, residual = 0.0016073271640315529
LinearSolver iteration 50, residual = 0.0010691425353595155
LinearSolver iteration 51, residual = 0.0010456470058024623
LinearSolver iteration 52, residual = 0.000595759153615959
LinearSolver iteration 53, residual = 0.0005939215989897622
LinearSolver iteration 54, residual = 0.00040290386616996244
LinearSolver iteration 55, residual = 0.00040066380782038003
LinearSolver iteration 56, residual = 0.0002543005761819086
LinearSolver iteration 57, residual = 0.0002527535623158061
LinearSolver iteration 58, residual = 0.00013526924429052613
LinearSolver iteration 59, residual = 0.00013376572240051968
LinearSolver iteration 60, residual = 5.1163604302161335e-05
LinearSolver iteration 61, residual = 5.0942439178644265e-05
LinearSolver iteration 62, residual = 2.443219452273615e-05
LinearSolver iteration 63, residual = 2.4186932007682973e-05
LinearSolver iteration 64, residual = 1.113132543105798e-05
LinearSolver iteration 65, residual = 1.1128361647134322e-05
LinearSolver iteration 66, residual = 4.244148442199896e-06
LinearSolver iteration 67, residual = 4.2413246789289464e-06
LinearSolver iteration 68, residual = 1.8252703488600543e-06
LinearSolver iteration 69, residual = 1.8248987304295922e-06
LinearSolver iteration 70, residual = 6.602943285291293e-07
LinearSolver iteration 71, residual = 6.602098965276247e-07
LinearSolver iteration 72, residual = 2.049931114192133e-07

basevector

Draw (gfu)

BaseWebGuiScene