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.320513884458098
MinRes iteration 3, residual = 1.9304918812163108
MinRes iteration 4, residual = 1.53396346307603
MinRes iteration 5, residual = 1.4866157469493735
MinRes iteration 6, residual = 1.277428164354539
MinRes iteration 7, residual = 1.2336050111420567
MinRes iteration 8, residual = 1.076716780243328
MinRes iteration 9, residual = 0.9764208183923647
MinRes iteration 10, residual = 0.8124102954934335
MinRes iteration 11, residual = 0.7630326288072056
MinRes iteration 12, residual = 0.6710519572747474
MinRes iteration 13, residual = 0.6468552815219445
MinRes iteration 14, residual = 0.5965412205029479
MinRes iteration 15, residual = 0.5729173657364667
MinRes iteration 16, residual = 0.5322237098138309
MinRes iteration 17, residual = 0.5244109437623257
MinRes iteration 18, residual = 0.4912002378193336
MinRes iteration 19, residual = 0.481630255599115
MinRes iteration 20, residual = 0.45495924807780896
MinRes iteration 21, residual = 0.4353886755547701
MinRes iteration 22, residual = 0.3922054493950942
MinRes iteration 23, residual = 0.37027788321903615
MinRes iteration 24, residual = 0.3011003348307442
MinRes iteration 25, residual = 0.2892067537855299
MinRes iteration 26, residual = 0.17308383727237375
MinRes iteration 27, residual = 0.16928179934916993
MinRes iteration 28, residual = 0.08657441731942168
MinRes iteration 29, residual = 0.08562898863022088
MinRes iteration 30, residual = 0.046854452342590065
MinRes iteration 31, residual = 0.04660020674568093
MinRes iteration 32, residual = 0.023441789137758386
MinRes iteration 33, residual = 0.02011841137302951
MinRes iteration 34, residual = 0.013048209565300133
MinRes iteration 35, residual = 0.012480171506024275
MinRes iteration 36, residual = 0.00913530597535238
MinRes iteration 37, residual = 0.008779885649184277
MinRes iteration 38, residual = 0.006460732446738416
MinRes iteration 39, residual = 0.006058729455792406
MinRes iteration 40, residual = 0.004153069889944102
MinRes iteration 41, residual = 0.003977244990868113
MinRes iteration 42, residual = 0.0025498351992951345
MinRes iteration 43, residual = 0.0023629754253481554
MinRes iteration 44, residual = 0.0014897024969477391
MinRes iteration 45, residual = 0.001438008170865018
MinRes iteration 46, residual = 0.0009456247234570568
MinRes iteration 47, residual = 0.000904225507869915
MinRes iteration 48, residual = 0.0006050248304509172
MinRes iteration 49, residual = 0.0005650216216937558
MinRes iteration 50, residual = 0.0003967204344211373
MinRes iteration 51, residual = 0.00037103618432571516
MinRes iteration 52, residual = 0.00024447080287521607
MinRes iteration 53, residual = 0.00023962944985242676
MinRes iteration 54, residual = 0.00013810579331522535
MinRes iteration 55, residual = 0.0001260702992063813
MinRes iteration 56, residual = 6.197735441207063e-05
MinRes iteration 57, residual = 5.419189162628389e-05
MinRes iteration 58, residual = 2.681788263678422e-05
MinRes iteration 59, residual = 2.3018259029955334e-05
MinRes iteration 60, residual = 1.0743945032634237e-05
MinRes iteration 61, residual = 1.0198694986634946e-05
MinRes iteration 62, residual = 4.66505716373943e-06
MinRes iteration 63, residual = 4.369820590534996e-06
MinRes iteration 64, residual = 1.877013071575374e-06
MinRes iteration 65, residual = 1.5246977870815582e-06
MinRes iteration 66, residual = 6.167769748720154e-07
MinRes iteration 67, residual = 5.436025690541738e-07
MinRes iteration 68, residual = 2.5093238287018675e-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.000179411319214
MinRes iteration 2, residual = 3.0834350054463613
MinRes iteration 3, residual = 2.116486149620493
MinRes iteration 4, residual = 1.8366080623996361
MinRes iteration 5, residual = 1.5895645511165304
MinRes iteration 6, residual = 1.544171088592255
MinRes iteration 7, residual = 1.4644985020244865
MinRes iteration 8, residual = 1.3521655569577125
MinRes iteration 9, residual = 1.3361032449540493
MinRes iteration 10, residual = 1.259562544285705
MinRes iteration 11, residual = 1.149502417016321
MinRes iteration 12, residual = 1.1491221325036691
MinRes iteration 13, residual = 1.0583582852278846
MinRes iteration 14, residual = 0.9773598464388681
MinRes iteration 15, residual = 0.9756045310923609
MinRes iteration 16, residual = 0.8670011646769751
MinRes iteration 17, residual = 0.7989406664024937
MinRes iteration 18, residual = 0.7843579390403238
MinRes iteration 19, residual = 0.6959411722935307
MinRes iteration 20, residual = 0.6804255085861398
MinRes iteration 21, residual = 0.6616381813443337
MinRes iteration 22, residual = 0.6248498919661581
MinRes iteration 23, residual = 0.6191429972376714
MinRes iteration 24, residual = 0.606798352966114
MinRes iteration 25, residual = 0.5794521829418521
MinRes iteration 26, residual = 0.5762377937624359
MinRes iteration 27, residual = 0.5583795182268426
MinRes iteration 28, residual = 0.5344331575565248
MinRes iteration 29, residual = 0.5342817943561792
MinRes iteration 30, residual = 0.5139643507350373
MinRes iteration 31, residual = 0.4969832846104731
MinRes iteration 32, residual = 0.4959591001669671
MinRes iteration 33, residual = 0.4768229915072394
MinRes iteration 34, residual = 0.4730163146750781
MinRes iteration 35, residual = 0.46161527016686277
MinRes iteration 36, residual = 0.4418922148215026
MinRes iteration 37, residual = 0.4418467069187686
MinRes iteration 38, residual = 0.42122052789830594
MinRes iteration 39, residual = 0.3993819042578253
MinRes iteration 40, residual = 0.3983788970407623
MinRes iteration 41, residual = 0.36309821165343165
MinRes iteration 42, residual = 0.34841248012378684
MinRes iteration 43, residual = 0.3264000812843558
MinRes iteration 44, residual = 0.2764309783815125
MinRes iteration 45, residual = 0.27642588506640214
MinRes iteration 46, residual = 0.22603268049146882
MinRes iteration 47, residual = 0.19615345110067275
MinRes iteration 48, residual = 0.18637317691621955
MinRes iteration 49, residual = 0.13823609788169272
MinRes iteration 50, residual = 0.13250126874036644
MinRes iteration 51, residual = 0.11243751436674804
MinRes iteration 52, residual = 0.08207853067976241
MinRes iteration 53, residual = 0.08173368881246544
MinRes iteration 54, residual = 0.060397357638124885
MinRes iteration 55, residual = 0.04495603845786887
MinRes iteration 56, residual = 0.043333358437040795
MinRes iteration 57, residual = 0.02691243354566012
MinRes iteration 58, residual = 0.023795887364581162
MinRes iteration 59, residual = 0.019506048626320948
MinRes iteration 60, residual = 0.013335974509605676
MinRes iteration 61, residual = 0.013253249776557669
MinRes iteration 62, residual = 0.009730686991315678
MinRes iteration 63, residual = 0.007615282174490367
MinRes iteration 64, residual = 0.007201992059253478
MinRes iteration 65, residual = 0.005138225367463208
MinRes iteration 66, residual = 0.005061934778905146
MinRes iteration 67, residual = 0.004238352793212358
MinRes iteration 68, residual = 0.0037078024627893555
MinRes iteration 69, residual = 0.003631483470068156
MinRes iteration 70, residual = 0.002832986896820194
MinRes iteration 71, residual = 0.0026378867956985993
MinRes iteration 72, residual = 0.0021883498902144053
MinRes iteration 73, residual = 0.0015223551846431386
MinRes iteration 74, residual = 0.0015223549158595291
MinRes iteration 75, residual = 0.0010546152670528584
MinRes iteration 76, residual = 0.0008252627892356858
MinRes iteration 77, residual = 0.0007549431451886332
MinRes iteration 78, residual = 0.0004634434695629877
MinRes iteration 79, residual = 0.0004556882337962502
MinRes iteration 80, residual = 0.0003038135607091898
MinRes iteration 81, residual = 0.00023835015421515748
MinRes iteration 82, residual = 0.00022166757617985374
MinRes iteration 83, residual = 0.00014830085987903085
MinRes iteration 84, residual = 0.00013126957452359918
MinRes iteration 85, residual = 0.00011184373407674518
MinRes iteration 86, residual = 6.946721782359395e-05
MinRes iteration 87, residual = 6.933932318997943e-05
MinRes iteration 88, residual = 4.494172974105446e-05
MinRes iteration 89, residual = 3.667501018554229e-05
MinRes iteration 90, residual = 3.156978437156566e-05
MinRes iteration 91, residual = 2.018745327196689e-05
MinRes iteration 92, residual = 1.964381468679682e-05
MinRes iteration 93, residual = 1.3523344178128678e-05
MinRes iteration 94, residual = 9.02921363335329e-06
MinRes iteration 95, residual = 8.616089336756953e-06
MinRes iteration 96, residual = 4.8851296082357605e-06
MinRes iteration 97, residual = 4.232314121545243e-06
MinRes iteration 98, residual = 3.2553482352209855e-06
MinRes iteration 99, residual = 1.929110390821917e-06
MinRes iteration 100, residual = 1.9249234530237667e-06
MinRes iteration 101, residual = 1.1575105286461877e-06
MinRes iteration 102, residual = 8.292932239162003e-07
MinRes iteration 103, residual = 7.279182383336709e-07
MinRes iteration 104, residual = 3.9761720422043324e-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.509846326863497
MinRes iteration 3, residual = 2.471375345209124
MinRes iteration 4, residual = 1.600041654102239
MinRes iteration 5, residual = 1.5999941924014434
MinRes iteration 6, residual = 1.331848127071812
MinRes iteration 7, residual = 1.2842023741344206
MinRes iteration 8, residual = 1.097267216042194
MinRes iteration 9, residual = 1.0536946299067136
MinRes iteration 10, residual = 0.8591400365732288
MinRes iteration 11, residual = 0.8019583070152935
MinRes iteration 12, residual = 0.6847419608924613
MinRes iteration 13, residual = 0.6763718338307138
MinRes iteration 14, residual = 0.6146066572025687
MinRes iteration 15, residual = 0.6023824596733502
MinRes iteration 16, residual = 0.5555606514054643
MinRes iteration 17, residual = 0.5430937217136549
MinRes iteration 18, residual = 0.5094093630083807
MinRes iteration 19, residual = 0.4995804424659091
MinRes iteration 20, residual = 0.4687233539171243
MinRes iteration 21, residual = 0.4561019247306363
MinRes iteration 22, residual = 0.42459537484596055
MinRes iteration 23, residual = 0.40438140563846586
MinRes iteration 24, residual = 0.3455546509234682
MinRes iteration 25, residual = 0.3189018314414594
MinRes iteration 26, residual = 0.2405341199646593
MinRes iteration 27, residual = 0.23388758308238228
MinRes iteration 28, residual = 0.13194535898020662
MinRes iteration 29, residual = 0.12123099213408245
MinRes iteration 30, residual = 0.06382557163323076
MinRes iteration 31, residual = 0.06303180309279424
MinRes iteration 32, residual = 0.030849945699984272
MinRes iteration 33, residual = 0.02934698171425437
MinRes iteration 34, residual = 0.015627485603871515
MinRes iteration 35, residual = 0.014617906806417256
MinRes iteration 36, residual = 0.008292871905733525
MinRes iteration 37, residual = 0.008150714306779407
MinRes iteration 38, residual = 0.005647994616945016
MinRes iteration 39, residual = 0.0056207663196512535
MinRes iteration 40, residual = 0.0038042890072358775
MinRes iteration 41, residual = 0.0038014749200140424
MinRes iteration 42, residual = 0.002091167783896991
MinRes iteration 43, residual = 0.002068155419536774
MinRes iteration 44, residual = 0.0010343859929532657
MinRes iteration 45, residual = 0.000945083731621424
MinRes iteration 46, residual = 0.0005040197865836424
MinRes iteration 47, residual = 0.0004947091691000631
MinRes iteration 48, residual = 0.00026577591276729837
MinRes iteration 49, residual = 0.00026255585601250815
MinRes iteration 50, residual = 0.000151696416517861
MinRes iteration 51, residual = 0.00015046497706185623
MinRes iteration 52, residual = 9.099101333534007e-05
MinRes iteration 53, residual = 9.064467690440782e-05
MinRes iteration 54, residual = 5.549445645770583e-05
MinRes iteration 55, residual = 5.44255529218227e-05
MinRes iteration 56, residual = 3.374035784884605e-05
MinRes iteration 57, residual = 3.3165537146835303e-05
MinRes iteration 58, residual = 1.8747435764210996e-05
MinRes iteration 59, residual = 1.8611134264364636e-05
MinRes iteration 60, residual = 9.84861610201104e-06
MinRes iteration 61, residual = 9.695503799562037e-06
MinRes iteration 62, residual = 5.197401872670218e-06
MinRes iteration 63, residual = 5.1835344118463265e-06
MinRes iteration 64, residual = 2.3650909689663154e-06
MinRes iteration 65, residual = 2.301651524107982e-06
MinRes iteration 66, residual = 8.885736555412273e-07

Draw (gfu)
Draw (gfp);