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.occ import *

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
Draw (shape);

geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
mesh.Curve(3)
Draw (mesh);

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

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

(u,p),(v,q) = X.TnT()

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

Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

gf = GridFunction(X)
gfu, gfp = gf.components

uin = CF ( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
Draw(gfu, mesh, min=0, max=2)
SetVisualization(max=2)

Solve equation:

res = -a.mat * gf.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
gf.vec.data += inv * res
Draw(gfu, mesh);

Testing different velocity-pressure pairs

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

def SolveStokes(X):
    (u,p),(v,q) = X.TnT()

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

    gf = GridFunction(X)
    gfu, gfp = gf.components

    uin = CF ( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
    gfu.Set(uin, definedon=mesh.Boundaries("inlet"))

    res = -a.mat * gf.vec
    inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
    gf.vec.data += inv * res

    Draw(gfu)

    return gfu

Higher order Taylor-Hood elements:

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

gfu = SolveStokes(X)

With discontinuous pressure elements P2-P1 is unstable:

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

gfu = SolveStokes(X)

V.ndof = 3368 , Q.ndof = 2364

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Cell In[9], line 6
      2 Q = L2(mesh, order=1)
      3 print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
      4 X = V*Q
      5
----> 6 gfu = SolveStokes(X)

Cell In[7], line 14, in SolveStokes(X)
     10     uin = CF ( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
     11     gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
     12
     13     res = -a.mat * gf.vec
---> 14     inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
     15     gf.vec.data += inv * res
     16
     17     Draw(gfu)

NgException: UmfpackInverse: Numeric factorization failed.

However, P2-P1 elements work on sub-divided simplicial meshes using Alfeld splits:

mesh2 = Mesh(geo.GenerateMesh(maxh=0.05)).Curve(3)
mesh2.SplitElements_Alfeld()
V = VectorH1(mesh2, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh2, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = V*Q
gfu = SolveStokes(X)

V.ndof = 9672 , Q.ndof = 7092

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

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

gfu = SolveStokes(X)

4944
V.ndof = 4944 , Q.ndof = 2364

the mini element:

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

gfu = SolveStokes(X)

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(InnerProduct(Grad(u),Grad(v))*dx).Assemble()
b = BilinearForm(div(u)*q*dx).Assemble()

Needed as preconditioner for the pressure:

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

Two right hand sides for the two spaces:

f = LinearForm(V).Assemble()
g = LinearForm(Q).Assemble();

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, printrates=True, initialize=False);

MinRes iteration 1, residual = 4.608947775183305
MinRes iteration 2, residual = 2.320513884458099
MinRes iteration 3, residual = 1.9304918812163114
MinRes iteration 4, residual = 1.5339634630760306
MinRes iteration 5, residual = 1.486615746949374
MinRes iteration 6, residual = 1.2774281643545398
MinRes iteration 7, residual = 1.2336050111420573
MinRes iteration 8, residual = 1.076716780243329
MinRes iteration 9, residual = 0.9764208183923657
MinRes iteration 10, residual = 0.8124102954934346
MinRes iteration 11, residual = 0.7630326288072061
MinRes iteration 12, residual = 0.6710519572747474
MinRes iteration 13, residual = 0.6468552815219438
MinRes iteration 14, residual = 0.5965412205029463
MinRes iteration 15, residual = 0.5729173657364642
MinRes iteration 16, residual = 0.532223709813827
MinRes iteration 17, residual = 0.5244109437623217
MinRes iteration 18, residual = 0.4912002378193273
MinRes iteration 19, residual = 0.48163025559910866
MinRes iteration 20, residual = 0.45495924807780086
MinRes iteration 21, residual = 0.4353886755547611
MinRes iteration 22, residual = 0.39220544939508234
MinRes iteration 23, residual = 0.37027788321902405
MinRes iteration 24, residual = 0.3011003348307297
MinRes iteration 25, residual = 0.2892067537855159
MinRes iteration 26, residual = 0.1730838372723632
MinRes iteration 27, residual = 0.16928179934915924
MinRes iteration 28, residual = 0.08657441731941727
MinRes iteration 29, residual = 0.08562898863021562
MinRes iteration 30, residual = 0.04685445234258631
MinRes iteration 31, residual = 0.04660020674567729
MinRes iteration 32, residual = 0.0234417891377562
MinRes iteration 33, residual = 0.020118411373027653
MinRes iteration 34, residual = 0.01304820956529891
MinRes iteration 35, residual = 0.012480171506023078
MinRes iteration 36, residual = 0.009135305975351506
MinRes iteration 37, residual = 0.008779885649183467
MinRes iteration 38, residual = 0.006460732446737782
MinRes iteration 39, residual = 0.006058729455791862
MinRes iteration 40, residual = 0.004153069889943694
MinRes iteration 41, residual = 0.0039772449908677555
MinRes iteration 42, residual = 0.002549835199294869
MinRes iteration 43, residual = 0.002362975425347911
MinRes iteration 44, residual = 0.0014897024969475637
MinRes iteration 45, residual = 0.0014380081708648418
MinRes iteration 46, residual = 0.0009456247234569397
MinRes iteration 47, residual = 0.000904225507869797
MinRes iteration 48, residual = 0.0006050248304508396
MinRes iteration 49, residual = 0.0005650216216936848
MinRes iteration 50, residual = 0.00039672043442108416
MinRes iteration 51, residual = 0.0003710361843256665
MinRes iteration 52, residual = 0.0002444708028751828
MinRes iteration 53, residual = 0.00023962944985239396
MinRes iteration 54, residual = 0.0001381057933152057
MinRes iteration 55, residual = 0.000126070299206362
MinRes iteration 56, residual = 6.197735441206218e-05
MinRes iteration 57, residual = 5.419189162627566e-05
MinRes iteration 58, residual = 2.6817882636780522e-05
MinRes iteration 59, residual = 2.3018259029951987e-05
MinRes iteration 60, residual = 1.074394503263267e-05
MinRes iteration 61, residual = 1.0198694986633389e-05
MinRes iteration 62, residual = 4.665057163738775e-06
MinRes iteration 63, residual = 4.369820590534393e-06
MinRes iteration 64, residual = 1.8770130715750827e-06
MinRes iteration 65, residual = 1.5246977870812935e-06
MinRes iteration 66, residual = 6.167769748719126e-07
MinRes iteration 67, residual = 5.436025690540513e-07
MinRes iteration 68, residual = 2.5093238287015197e-07

Draw (gfu);

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(p*q*dx).Assemble()

f = LinearForm(V)
f.Assemble()

g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")
uin = CF( (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, printrates=True, initialize=False, maxsteps=200);

MinRes iteration 1, residual = 7.0001794113192135
MinRes iteration 2, residual = 3.0834350054463604
MinRes iteration 3, residual = 2.11648614962049
MinRes iteration 4, residual = 1.8366080623996337
MinRes iteration 5, residual = 1.5895645511165277
MinRes iteration 6, residual = 1.5441710885922524
MinRes iteration 7, residual = 1.464498502024484
MinRes iteration 8, residual = 1.3521655569577107
MinRes iteration 9, residual = 1.3361032449540478
MinRes iteration 10, residual = 1.259562544285704
MinRes iteration 11, residual = 1.1495024170163202
MinRes iteration 12, residual = 1.1491221325036682
MinRes iteration 13, residual = 1.058358285227884
MinRes iteration 14, residual = 0.977359846438868
MinRes iteration 15, residual = 0.9756045310923608
MinRes iteration 16, residual = 0.8670011646769761
MinRes iteration 17, residual = 0.7989406664024945
MinRes iteration 18, residual = 0.7843579390403249
MinRes iteration 19, residual = 0.6959411722935319
MinRes iteration 20, residual = 0.6804255085861413
MinRes iteration 21, residual = 0.6616381813443353
MinRes iteration 22, residual = 0.6248498919661603
MinRes iteration 23, residual = 0.6191429972376736
MinRes iteration 24, residual = 0.6067983529661163
MinRes iteration 25, residual = 0.5794521829418546
MinRes iteration 26, residual = 0.5762377937624386
MinRes iteration 27, residual = 0.5583795182268455
MinRes iteration 28, residual = 0.5344331575565282
MinRes iteration 29, residual = 0.5342817943561826
MinRes iteration 30, residual = 0.5139643507350413
MinRes iteration 31, residual = 0.4969832846104772
MinRes iteration 32, residual = 0.4959591001669713
MinRes iteration 33, residual = 0.47682299150724383
MinRes iteration 34, residual = 0.4730163146750826
MinRes iteration 35, residual = 0.4616152701668676
MinRes iteration 36, residual = 0.4418922148215079
MinRes iteration 37, residual = 0.4418467069187739
MinRes iteration 38, residual = 0.4212205278983119
MinRes iteration 39, residual = 0.3993819042578315
MinRes iteration 40, residual = 0.39837889704076856
MinRes iteration 41, residual = 0.36309821165343836
MinRes iteration 42, residual = 0.34841248012379306
MinRes iteration 43, residual = 0.32640008128436315
MinRes iteration 44, residual = 0.27643097838152053
MinRes iteration 45, residual = 0.27642588506641025
MinRes iteration 46, residual = 0.22603268049149294
MinRes iteration 47, residual = 0.1961534511006955
MinRes iteration 48, residual = 0.18637317691636587
MinRes iteration 49, residual = 0.1382360978821064
MinRes iteration 50, residual = 0.13250126874040338
MinRes iteration 51, residual = 0.11243751436736768
MinRes iteration 52, residual = 0.08207853067984679
MinRes iteration 53, residual = 0.08173368881251641
MinRes iteration 54, residual = 0.060397357638157526
MinRes iteration 55, residual = 0.0449560384578713
MinRes iteration 56, residual = 0.04333335843704772
MinRes iteration 57, residual = 0.0269124335456661
MinRes iteration 58, residual = 0.023795887364581908
MinRes iteration 59, residual = 0.019506048626337615
MinRes iteration 60, residual = 0.013335974509627513
MinRes iteration 61, residual = 0.013253249776568596
MinRes iteration 62, residual = 0.009730686991462468
MinRes iteration 63, residual = 0.007615282174583801
MinRes iteration 64, residual = 0.0072019920598046785
MinRes iteration 65, residual = 0.005138225368955788
MinRes iteration 66, residual = 0.00506193477924885
MinRes iteration 67, residual = 0.004238352800437373
MinRes iteration 68, residual = 0.003707802469130561
MinRes iteration 69, residual = 0.003631483497045454
MinRes iteration 70, residual = 0.0028329870828225955
MinRes iteration 71, residual = 0.0026378867968233585
MinRes iteration 72, residual = 0.0021883513713352866
MinRes iteration 73, residual = 0.0015223581700689827
MinRes iteration 74, residual = 0.001522357890988066
MinRes iteration 75, residual = 0.0010546477312166817
MinRes iteration 76, residual = 0.0008252820541971611
MinRes iteration 77, residual = 0.0007551068353856259
MinRes iteration 78, residual = 0.00046390346194408796
MinRes iteration 79, residual = 0.00045584458538497195
MinRes iteration 80, residual = 0.0003063016314149405
MinRes iteration 81, residual = 0.00023899387877681657
MinRes iteration 82, residual = 0.00022613886276693984
MinRes iteration 83, residual = 0.00015293656996993413
MinRes iteration 84, residual = 0.00013132846208975142
MinRes iteration 85, residual = 0.00011320892153436389
MinRes iteration 86, residual = 6.956482874531344e-05
MinRes iteration 87, residual = 6.941163033838303e-05
MinRes iteration 88, residual = 4.496026617909861e-05
MinRes iteration 89, residual = 3.6675168372936237e-05
MinRes iteration 90, residual = 3.1571575141521166e-05
MinRes iteration 91, residual = 2.0187642755748624e-05
MinRes iteration 92, residual = 1.9643874259209198e-05
MinRes iteration 93, residual = 1.3523379440721044e-05
MinRes iteration 94, residual = 9.029213949021859e-06
MinRes iteration 95, residual = 8.61609182855131e-06
MinRes iteration 96, residual = 4.885130019371125e-06
MinRes iteration 97, residual = 4.2323141429097565e-06
MinRes iteration 98, residual = 3.255348458878199e-06
MinRes iteration 99, residual = 1.9291108192019626e-06
MinRes iteration 100, residual = 1.9249237521671566e-06
MinRes iteration 101, residual = 1.1575158303848566e-06
MinRes iteration 102, residual = 8.292964023922225e-07
MinRes iteration 103, residual = 7.279486066484861e-07
MinRes iteration 104, residual = 3.9771388338026166e-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(p*q*dx).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.Assemble()

g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")
uin = CF( (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=True, initialize=False, maxsteps=200);

MinRes iteration 1, residual = 10.072431117077517
MinRes iteration 2, residual = 3.5098463268634976
MinRes iteration 3, residual = 2.471375345209127
MinRes iteration 4, residual = 1.600041654102241
MinRes iteration 5, residual = 1.5999941924014456
MinRes iteration 6, residual = 1.3318481270718132
MinRes iteration 7, residual = 1.2842023741344215
MinRes iteration 8, residual = 1.097267216042194
MinRes iteration 9, residual = 1.0536946299067131
MinRes iteration 10, residual = 0.8591400365732266
MinRes iteration 11, residual = 0.801958307015291
MinRes iteration 12, residual = 0.684741960892458
MinRes iteration 13, residual = 0.6763718338307104
MinRes iteration 14, residual = 0.6146066572025645
MinRes iteration 15, residual = 0.6023824596733459
MinRes iteration 16, residual = 0.555560651405459
MinRes iteration 17, residual = 0.5430937217136498
MinRes iteration 18, residual = 0.5094093630083749
MinRes iteration 19, residual = 0.4995804424659035
MinRes iteration 20, residual = 0.4687233539171181
MinRes iteration 21, residual = 0.45610192473063055
MinRes iteration 22, residual = 0.4245953748459544
MinRes iteration 23, residual = 0.40438140563846026
MinRes iteration 24, residual = 0.34555465092346194
MinRes iteration 25, residual = 0.3189018314414527
MinRes iteration 26, residual = 0.24053411996465315
MinRes iteration 27, residual = 0.2338875830823767
MinRes iteration 28, residual = 0.13194535898019727
MinRes iteration 29, residual = 0.12123099213407908
MinRes iteration 30, residual = 0.06382557163322844
MinRes iteration 31, residual = 0.06303180309279277
MinRes iteration 32, residual = 0.030849945699983068
MinRes iteration 33, residual = 0.02934698171425281
MinRes iteration 34, residual = 0.01562748560387062
MinRes iteration 35, residual = 0.014617906806416154
MinRes iteration 36, residual = 0.008292871905733048
MinRes iteration 37, residual = 0.008150714306778939
MinRes iteration 38, residual = 0.005647994616944674
MinRes iteration 39, residual = 0.005620766319650907
MinRes iteration 40, residual = 0.003804289007235603
MinRes iteration 41, residual = 0.003801474920013776
MinRes iteration 42, residual = 0.0020911677838968555
MinRes iteration 43, residual = 0.0020681554195366297
MinRes iteration 44, residual = 0.0010343859929531985
MinRes iteration 45, residual = 0.0009450837316213513
MinRes iteration 46, residual = 0.0005040197865836086
MinRes iteration 47, residual = 0.0004947091691000279
MinRes iteration 48, residual = 0.00026577591276728086
MinRes iteration 49, residual = 0.0002625558560124905
MinRes iteration 50, residual = 0.0001516964165178512
MinRes iteration 51, residual = 0.0001504649770618458
MinRes iteration 52, residual = 9.099101333533473e-05
MinRes iteration 53, residual = 9.064467690440229e-05
MinRes iteration 54, residual = 5.549445645770321e-05
MinRes iteration 55, residual = 5.4425552921819985e-05
MinRes iteration 56, residual = 3.3740357848844794e-05
MinRes iteration 57, residual = 3.316553714683386e-05
MinRes iteration 58, residual = 1.874743576421012e-05
MinRes iteration 59, residual = 1.8611134264363718e-05
MinRes iteration 60, residual = 9.848616102011318e-06
MinRes iteration 61, residual = 9.69550379956231e-06
MinRes iteration 62, residual = 5.197401872672237e-06
MinRes iteration 63, residual = 5.183534411848338e-06
MinRes iteration 64, residual = 2.3650909689684952e-06
MinRes iteration 65, residual = 2.301651524109962e-06
MinRes iteration 66, residual = 8.885736555422786e-07

Draw (gfu)
Draw (gfp);