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 0x7fa12e7b9f70>

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)
/tmp/ipykernel_31924/1189186610.py in <module>
      4 X = V*V*Q
      5
----> 6 gfu = SolveStokes(X)

/tmp/ipykernel_31924/811579872.py 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 0x7fa12e9d0930>

Needed as preconditioner for the pressure:

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

<ngsolve.comp.BilinearForm at 0x7fa12e7dddb0>

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 0x7fa12e7db930>

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

A stabilised lowest order formulation:

mesh = Mesh( geo.GenerateMesh(maxh=0.02))

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

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

a = BilinearForm(InnerProduct(Grad(u),Grad(v))*dx).Assemble()
b = BilinearForm(div(u)*q*dx).Assemble()
h = specialcf.mesh_size
c = BilinearForm(-0.1*h*h*grad(p)*grad(q)*dx).Assemble()

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

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

g = LinearForm(Q)
g.Assemble()

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"))

K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, c.mat] ] )
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, maxsteps=200);

LinearSolver iteration 1, residual = 7.093524574559076
LinearSolver iteration 2, residual = 3.123301846060359
LinearSolver iteration 3, residual = 2.1396216166041753
LinearSolver iteration 4, residual = 1.914219462615039
LinearSolver iteration 5, residual = 1.6116196031976782
LinearSolver iteration 6, residual = 1.5489603152856475
LinearSolver iteration 7, residual = 1.4848591185063968
LinearSolver iteration 8, residual = 1.3800580493063581
LinearSolver iteration 9, residual = 1.3386235485379536
LinearSolver iteration 10, residual = 1.284152849911857
LinearSolver iteration 11, residual = 1.171400089134812
LinearSolver iteration 12, residual = 1.1640324250820986
LinearSolver iteration 13, residual = 1.0888266346366466
LinearSolver iteration 14, residual = 1.0126085439860515
LinearSolver iteration 15, residual = 1.0008115366910437
LinearSolver iteration 16, residual = 0.8910955738823616
LinearSolver iteration 17, residual = 0.8729880985957466
LinearSolver iteration 18, residual = 0.7981674440668808
LinearSolver iteration 19, residual = 0.7244171610865654
LinearSolver iteration 20, residual = 0.7233058879256861
LinearSolver iteration 21, residual = 0.6687530797250423
LinearSolver iteration 22, residual = 0.6526056449597718
LinearSolver iteration 23, residual = 0.6436431389122804
LinearSolver iteration 24, residual = 0.6159761414220264
LinearSolver iteration 25, residual = 0.6115222527483778
LinearSolver iteration 26, residual = 0.5991667862423437
LinearSolver iteration 27, residual = 0.5736793196432238
LinearSolver iteration 28, residual = 0.5714962016331891
LinearSolver iteration 29, residual = 0.5542631844122419
LinearSolver iteration 30, residual = 0.532055024445272
LinearSolver iteration 31, residual = 0.531921457390652
LinearSolver iteration 32, residual = 0.510729746762853
LinearSolver iteration 33, residual = 0.5019586208536521
LinearSolver iteration 34, residual = 0.49560770088990036
LinearSolver iteration 35, residual = 0.4750072960858663
LinearSolver iteration 36, residual = 0.4713090543262219
LinearSolver iteration 37, residual = 0.46146376709121284
LinearSolver iteration 38, residual = 0.4391832079345872
LinearSolver iteration 39, residual = 0.437368683586536
LinearSolver iteration 40, residual = 0.4211512919810195
LinearSolver iteration 41, residual = 0.39367115138774345
LinearSolver iteration 42, residual = 0.39347657129658536
LinearSolver iteration 43, residual = 0.36001180083278406
LinearSolver iteration 44, residual = 0.33039855672199453
LinearSolver iteration 45, residual = 0.3246152628977096
LinearSolver iteration 46, residual = 0.26709807045403816
LinearSolver iteration 47, residual = 0.25593084790674575
LinearSolver iteration 48, residual = 0.22792999737229336
LinearSolver iteration 49, residual = 0.18718749769798987
LinearSolver iteration 50, residual = 0.18567357579428906
LinearSolver iteration 51, residual = 0.1351835183435528
LinearSolver iteration 52, residual = 0.12299160854248277
LinearSolver iteration 53, residual = 0.10499773205792497
LinearSolver iteration 54, residual = 0.07415268727177486
LinearSolver iteration 55, residual = 0.07292909225916672
LinearSolver iteration 56, residual = 0.0534475562309993
LinearSolver iteration 57, residual = 0.03862119936371087
LinearSolver iteration 58, residual = 0.03670906203305091
LinearSolver iteration 59, residual = 0.023268422282318204
LinearSolver iteration 60, residual = 0.022768409341301463
LinearSolver iteration 61, residual = 0.0166865078791172
LinearSolver iteration 62, residual = 0.012517856611925551
LinearSolver iteration 63, residual = 0.012287318123687411
LinearSolver iteration 64, residual = 0.00835034915350182
LinearSolver iteration 65, residual = 0.007108799169281131
LinearSolver iteration 66, residual = 0.006255289192429071
LinearSolver iteration 67, residual = 0.004658268592092344
LinearSolver iteration 68, residual = 0.004614851095049302
LinearSolver iteration 69, residual = 0.003926067368038764
LinearSolver iteration 70, residual = 0.0033784110116021836
LinearSolver iteration 71, residual = 0.003339413114352952
LinearSolver iteration 72, residual = 0.0025256152680784155
LinearSolver iteration 73, residual = 0.0022577605196393885
LinearSolver iteration 74, residual = 0.0019389233894438346
LinearSolver iteration 75, residual = 0.0013265439371479171
LinearSolver iteration 76, residual = 0.0013047760787624745
LinearSolver iteration 77, residual = 0.0009821054137240845
LinearSolver iteration 78, residual = 0.0006879917754877875
LinearSolver iteration 79, residual = 0.0006866656904256515
LinearSolver iteration 80, residual = 0.0004627504518609887
LinearSolver iteration 81, residual = 0.0003366306396003939
LinearSolver iteration 82, residual = 0.0003247538495450309
LinearSolver iteration 83, residual = 0.00020472316694769625
LinearSolver iteration 84, residual = 0.00017302449351049687
LinearSolver iteration 85, residual = 0.00015055955198062793
LinearSolver iteration 86, residual = 9.503979483358908e-05
LinearSolver iteration 87, residual = 9.139967370081062e-05
LinearSolver iteration 88, residual = 6.205580470675907e-05
LinearSolver iteration 89, residual = 4.417482518300715e-05
LinearSolver iteration 90, residual = 4.208015931432418e-05
LinearSolver iteration 91, residual = 2.639415198204693e-05
LinearSolver iteration 92, residual = 2.3198211487033754e-05
LinearSolver iteration 93, residual = 2.05976446848021e-05
LinearSolver iteration 94, residual = 1.335750488850388e-05
LinearSolver iteration 95, residual = 1.178183346676212e-05
LinearSolver iteration 96, residual = 9.221492579476037e-06
LinearSolver iteration 97, residual = 5.430999939758362e-06
LinearSolver iteration 98, residual = 5.424479296528133e-06
LinearSolver iteration 99, residual = 3.2829817157318746e-06
LinearSolver iteration 100, residual = 2.3963406655547116e-06
LinearSolver iteration 101, residual = 2.1384549823450807e-06
LinearSolver iteration 102, residual = 1.1841165235066038e-06
LinearSolver iteration 103, residual = 1.1037166205473494e-06
LinearSolver iteration 104, residual = 7.415661381672827e-07
LinearSolver iteration 105, residual = 4.5723926977354625e-07

Draw (gfu)
Draw (gfp);

(P1-iso-P2)-P1

The mesh is uniformly refined once. The velocity space is \(P^1\) on the refined mesh, but the pressure is continuous \(P^1\) in the coarse space. We obtain this by representing the pressure in the range of the coarse-to-fine prolongation operator.

mesh = Mesh( geo.GenerateMesh(maxh=0.02))

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

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

mp = BilinearForm(Q)
mp += p*q*dx
mp.Assemble()
invmp = mp.mat.Inverse(inverse="sparsecholesky")

mesh.ngmesh.Refine()
V.Update()
Q.Update()

# prolongation operator from mesh-level 0 to level 1:
prol = Q.Prolongation().Operator(1)
invmp2 = prol @ invmp @ prol.T

a = BilinearForm(InnerProduct(Grad(u),Grad(v))*dx).Assemble()
b = BilinearForm(div(u)*q*dx).Assemble()

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

g = LinearForm(Q)
g.Assemble()

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"))

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

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

solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, \
                printrates='\r', initialize=False, maxsteps=200);

LinearSolver converged in 70 iterations to residual 6.656296451166797e-07

Draw (gfu)
Draw (gfp);