# 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:

[1]:

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

[2]:

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:

[3]:

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.

[4]:

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

a = BilinearForm(stokes).Assemble()


Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

[5]:

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:

[6]:

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:

[7]:

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

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:

[8]:

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:

[9]:

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 = 3424 , Q.ndof = 2406

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

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

NgException: UmfpackInverse: Numeric factorization failed.


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

[10]:

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 = 9840 , Q.ndof = 7218


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

[11]:

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)

5028
V.ndof = 5028 , Q.ndof = 2406


the mini element:

[12]:

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:

[13]:

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

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

b = BilinearForm(div(u)*q*dx).Assemble()


Needed as preconditioner for the pressure:

[14]:

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


Two right hand sides for the two spaces:

[15]:

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


Two GridFunctions for velocity and pressure:

[16]:

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.

[17]:

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

LinearSolver iteration 1, residual = 4.518937179332287
LinearSolver iteration 2, residual = 2.3020782066183405
LinearSolver iteration 3, residual = 1.9171806179628295
LinearSolver iteration 4, residual = 1.5318309921095468
LinearSolver iteration 5, residual = 1.4827182650518633
LinearSolver iteration 6, residual = 1.2760882423711855
LinearSolver iteration 7, residual = 1.2316328218317878
LinearSolver iteration 8, residual = 1.0758256973668074
LinearSolver iteration 9, residual = 0.9742346042735229
LinearSolver iteration 10, residual = 0.8116311768045531
LinearSolver iteration 11, residual = 0.7620621123363509
LinearSolver iteration 12, residual = 0.6711057838761307
LinearSolver iteration 13, residual = 0.6467734480369526
LinearSolver iteration 14, residual = 0.5963885550482779
LinearSolver iteration 15, residual = 0.5715760562984827
LinearSolver iteration 16, residual = 0.5314172872757961
LinearSolver iteration 17, residual = 0.523188572014509
LinearSolver iteration 18, residual = 0.48916660118256133
LinearSolver iteration 19, residual = 0.48026612300248067
LinearSolver iteration 20, residual = 0.45410105844800014
LinearSolver iteration 21, residual = 0.4339175906792935
LinearSolver iteration 22, residual = 0.38970115024248825
LinearSolver iteration 23, residual = 0.3665030310531088
LinearSolver iteration 24, residual = 0.29515807082938256
LinearSolver iteration 25, residual = 0.2840989870359845
LinearSolver iteration 26, residual = 0.17263506434502293
LinearSolver iteration 27, residual = 0.16533278871221124
LinearSolver iteration 28, residual = 0.09207002076023608
LinearSolver iteration 29, residual = 0.08806383521916354
LinearSolver iteration 30, residual = 0.045707773020785346
LinearSolver iteration 31, residual = 0.04558872526001847
LinearSolver iteration 32, residual = 0.022509806937679445
LinearSolver iteration 33, residual = 0.01911902981935344
LinearSolver iteration 34, residual = 0.012572671188669397
LinearSolver iteration 35, residual = 0.012232920688064624
LinearSolver iteration 36, residual = 0.008796244991863617
LinearSolver iteration 37, residual = 0.008540624451742073
LinearSolver iteration 38, residual = 0.006168751068326351
LinearSolver iteration 39, residual = 0.005859614741848614
LinearSolver iteration 40, residual = 0.003919721191004725
LinearSolver iteration 41, residual = 0.0037309231359978446
LinearSolver iteration 42, residual = 0.002409644869021594
LinearSolver iteration 43, residual = 0.002233246013776026
LinearSolver iteration 44, residual = 0.001440117277538793
LinearSolver iteration 45, residual = 0.0014169553919431872
LinearSolver iteration 46, residual = 0.0009214338668519829
LinearSolver iteration 47, residual = 0.0008998464054884886
LinearSolver iteration 48, residual = 0.0005762503566586648
LinearSolver iteration 49, residual = 0.0005472865054142066
LinearSolver iteration 50, residual = 0.0003758901019645419
LinearSolver iteration 51, residual = 0.0003524418202292021
LinearSolver iteration 52, residual = 0.00023429185631029594
LinearSolver iteration 53, residual = 0.000231577683020344
LinearSolver iteration 54, residual = 0.00013261675061463022
LinearSolver iteration 55, residual = 0.00012528252842354204
LinearSolver iteration 56, residual = 5.852579935273627e-05
LinearSolver iteration 57, residual = 5.220921278375006e-05
LinearSolver iteration 58, residual = 2.445981833622924e-05
LinearSolver iteration 59, residual = 2.1011457947562368e-05
LinearSolver iteration 60, residual = 9.894725207855162e-06
LinearSolver iteration 61, residual = 9.325347229145796e-06
LinearSolver iteration 62, residual = 4.39262521268901e-06
LinearSolver iteration 63, residual = 4.195388653320418e-06
LinearSolver iteration 64, residual = 1.7610904076427593e-06
LinearSolver iteration 65, residual = 1.487707702456541e-06
LinearSolver iteration 66, residual = 5.821123043092815e-07
LinearSolver iteration 67, residual = 5.153650968031574e-07
LinearSolver iteration 68, residual = 2.3534539744561125e-07

[18]:

Draw (gfu);


## A stabilised lowest order formulation:¶

[19]:

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

b = BilinearForm(div(u)*q*dx).Assemble()
h = specialcf.mesh_size

[20]:

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

LinearSolver iteration 1, residual = 7.103717015785737
LinearSolver iteration 2, residual = 3.1240529736545715
LinearSolver iteration 3, residual = 2.132679123744692
LinearSolver iteration 4, residual = 1.8523744559707935
LinearSolver iteration 5, residual = 1.5949503379856076
LinearSolver iteration 6, residual = 1.5480493182937696
LinearSolver iteration 7, residual = 1.4639703507160722
LinearSolver iteration 8, residual = 1.3511179026398876
LinearSolver iteration 9, residual = 1.3380334785486305
LinearSolver iteration 10, residual = 1.2598116080932027
LinearSolver iteration 11, residual = 1.1487182058160144
LinearSolver iteration 12, residual = 1.1487081802610593
LinearSolver iteration 13, residual = 1.0553572664456368
LinearSolver iteration 14, residual = 0.976593791636269
LinearSolver iteration 15, residual = 0.9718806996194856
LinearSolver iteration 16, residual = 0.8599661205136232
LinearSolver iteration 17, residual = 0.8005914699770541
LinearSolver iteration 18, residual = 0.7815317337264195
LinearSolver iteration 19, residual = 0.695091815685787
LinearSolver iteration 20, residual = 0.681097261054522
LinearSolver iteration 21, residual = 0.6601652993538863
LinearSolver iteration 22, residual = 0.6238986254835396
LinearSolver iteration 23, residual = 0.6195450703831556
LinearSolver iteration 24, residual = 0.6055276228835791
LinearSolver iteration 25, residual = 0.5784318059307156
LinearSolver iteration 26, residual = 0.5763221221836843
LinearSolver iteration 27, residual = 0.5571394201839083
LinearSolver iteration 28, residual = 0.5339897335990003
LinearSolver iteration 29, residual = 0.5339897224407189
LinearSolver iteration 30, residual = 0.5122111603632727
LinearSolver iteration 31, residual = 0.49709386368401864
LinearSolver iteration 32, residual = 0.49522266954245237
LinearSolver iteration 33, residual = 0.4757510640129083
LinearSolver iteration 34, residual = 0.4731826388330148
LinearSolver iteration 35, residual = 0.4601388105401166
LinearSolver iteration 36, residual = 0.4412886936788078
LinearSolver iteration 37, residual = 0.44120267904223176
LinearSolver iteration 38, residual = 0.4187428770452127
LinearSolver iteration 39, residual = 0.39925869154183435
LinearSolver iteration 40, residual = 0.3965165497660975
LinearSolver iteration 41, residual = 0.35958302389869107
LinearSolver iteration 42, residual = 0.34918705384313015
LinearSolver iteration 43, residual = 0.32081334730538613
LinearSolver iteration 44, residual = 0.27496623675114984
LinearSolver iteration 45, residual = 0.2739395154602775
LinearSolver iteration 46, residual = 0.22065410237349553
LinearSolver iteration 47, residual = 0.19723595462027937
LinearSolver iteration 48, residual = 0.18105301698811066
LinearSolver iteration 49, residual = 0.13692381837465906
LinearSolver iteration 50, residual = 0.13287558442797537
LinearSolver iteration 51, residual = 0.10926674275131437
LinearSolver iteration 52, residual = 0.08098202942409063
LinearSolver iteration 53, residual = 0.08094489870875758
LinearSolver iteration 54, residual = 0.056389394203058316
LinearSolver iteration 55, residual = 0.04514989779565494
LinearSolver iteration 56, residual = 0.040324817474304074
LinearSolver iteration 57, residual = 0.02515404088463173
LinearSolver iteration 58, residual = 0.023960521236286676
LinearSolver iteration 59, residual = 0.018012559701725786
LinearSolver iteration 60, residual = 0.01300583626679384
LinearSolver iteration 61, residual = 0.012973772374450137
LinearSolver iteration 62, residual = 0.009169395609592765
LinearSolver iteration 63, residual = 0.00763667146865482
LinearSolver iteration 64, residual = 0.0068485873124133475
LinearSolver iteration 65, residual = 0.004943988286513263
LinearSolver iteration 66, residual = 0.004911909351469885
LinearSolver iteration 67, residual = 0.004081919199534889
LinearSolver iteration 68, residual = 0.0035790619585567927
LinearSolver iteration 69, residual = 0.0034957389244193965
LinearSolver iteration 70, residual = 0.002702516402196035
LinearSolver iteration 71, residual = 0.0026133921922323797
LinearSolver iteration 72, residual = 0.002044199670776769
LinearSolver iteration 73, residual = 0.0014952781332316647
LinearSolver iteration 74, residual = 0.0014711480663878792
LinearSolver iteration 75, residual = 0.000991617872172705
LinearSolver iteration 76, residual = 0.0008714858820776782
LinearSolver iteration 77, residual = 0.0007056267834200474
LinearSolver iteration 78, residual = 0.0004773645624879543
LinearSolver iteration 79, residual = 0.00046424789573374227
LinearSolver iteration 80, residual = 0.00028094430610380993
LinearSolver iteration 81, residual = 0.0002548453786863786
LinearSolver iteration 82, residual = 0.00020256676252789458
LinearSolver iteration 83, residual = 0.00013267945882768027
LinearSolver iteration 84, residual = 0.00013194521868198242
LinearSolver iteration 85, residual = 8.593924383152689e-05
LinearSolver iteration 86, residual = 6.239009989375969e-05
LinearSolver iteration 87, residual = 5.9815516703039264e-05
LinearSolver iteration 88, residual = 3.849464222210089e-05
LinearSolver iteration 89, residual = 3.2605363484935703e-05
LinearSolver iteration 90, residual = 2.8322245608267862e-05
LinearSolver iteration 91, residual = 1.82690807881192e-05
LinearSolver iteration 92, residual = 1.8051591698283513e-05
LinearSolver iteration 93, residual = 1.1906749700372199e-05
LinearSolver iteration 94, residual = 8.918371402129329e-06
LinearSolver iteration 95, residual = 7.708397821999984e-06
LinearSolver iteration 96, residual = 4.360196894806923e-06
LinearSolver iteration 97, residual = 4.241381799516438e-06
LinearSolver iteration 98, residual = 2.799363495931647e-06
LinearSolver iteration 99, residual = 1.8482147124704675e-06
LinearSolver iteration 100, residual = 1.7777878936095626e-06
LinearSolver iteration 101, residual = 9.742813141308085e-07
LinearSolver iteration 102, residual = 8.494431817617426e-07
LinearSolver iteration 103, residual = 6.2071818814675e-07

[21]:

Draw (gfu);

[22]:

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.

[23]:

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

[24]:

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

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

LinearSolver iteration 1, residual = 10.21762787865306
LinearSolver iteration 2, residual = 3.5470429131109387
LinearSolver iteration 3, residual = 2.4903542579953064
LinearSolver iteration 4, residual = 1.6025072100687434
LinearSolver iteration 5, residual = 1.6024983010872687
LinearSolver iteration 6, residual = 1.3330581375717194
LinearSolver iteration 7, residual = 1.2854650025118093
LinearSolver iteration 8, residual = 1.0977943618957249
LinearSolver iteration 9, residual = 1.0548387055316293
LinearSolver iteration 10, residual = 0.8601293957401446
LinearSolver iteration 11, residual = 0.8024876442330512
LinearSolver iteration 12, residual = 0.6849389452406769
LinearSolver iteration 13, residual = 0.67661040458432
LinearSolver iteration 14, residual = 0.6147786664099295
LinearSolver iteration 15, residual = 0.602528668180872
LinearSolver iteration 16, residual = 0.5556519286013163
LinearSolver iteration 17, residual = 0.5433107795923491
LinearSolver iteration 18, residual = 0.5095539811129507
LinearSolver iteration 19, residual = 0.49972652714034743
LinearSolver iteration 20, residual = 0.4689365281768766
LinearSolver iteration 21, residual = 0.4562975098172982
LinearSolver iteration 22, residual = 0.4247374804729204
LinearSolver iteration 23, residual = 0.40470386462061436
LinearSolver iteration 24, residual = 0.346152345360998
LinearSolver iteration 25, residual = 0.319072865792631
LinearSolver iteration 26, residual = 0.24068756246455486
LinearSolver iteration 27, residual = 0.2340879104748355
LinearSolver iteration 28, residual = 0.1321028667368393
LinearSolver iteration 29, residual = 0.12149542221441378
LinearSolver iteration 30, residual = 0.06395535062746384
LinearSolver iteration 31, residual = 0.06313158950829824
LinearSolver iteration 32, residual = 0.03089930268379747
LinearSolver iteration 33, residual = 0.02941571373123378
LinearSolver iteration 34, residual = 0.01569096172127071
LinearSolver iteration 35, residual = 0.014673468623518305
LinearSolver iteration 36, residual = 0.008379750118117056
LinearSolver iteration 37, residual = 0.008238251565402985
LinearSolver iteration 38, residual = 0.005709345997468741
LinearSolver iteration 39, residual = 0.00568497523031095
LinearSolver iteration 40, residual = 0.003835868766848464
LinearSolver iteration 41, residual = 0.003832185229123109
LinearSolver iteration 42, residual = 0.002101456444473924
LinearSolver iteration 43, residual = 0.002079171810906028
LinearSolver iteration 44, residual = 0.0010349257957432318
LinearSolver iteration 45, residual = 0.0009452428556789443
LinearSolver iteration 46, residual = 0.0005018308254017286
LinearSolver iteration 47, residual = 0.0004925087878097115
LinearSolver iteration 48, residual = 0.00026249207049355356
LinearSolver iteration 49, residual = 0.00025892211488885697
LinearSolver iteration 50, residual = 0.0001490547702294028
LinearSolver iteration 51, residual = 0.0001476873822811908
LinearSolver iteration 52, residual = 8.918410229409894e-05
LinearSolver iteration 53, residual = 8.884336870666497e-05
LinearSolver iteration 54, residual = 5.443788069474884e-05
LinearSolver iteration 55, residual = 5.330171690491297e-05
LinearSolver iteration 56, residual = 3.318034823067684e-05
LinearSolver iteration 57, residual = 3.259166338176218e-05
LinearSolver iteration 58, residual = 1.8456984186068262e-05
LinearSolver iteration 59, residual = 1.8306187978097573e-05
LinearSolver iteration 60, residual = 9.729978738596325e-06
LinearSolver iteration 61, residual = 9.554842425624725e-06
LinearSolver iteration 62, residual = 5.161704605550565e-06
LinearSolver iteration 63, residual = 5.149417860117464e-06
LinearSolver iteration 64, residual = 2.365836031715049e-06
LinearSolver iteration 65, residual = 2.2957609093107117e-06
LinearSolver iteration 66, residual = 8.861914669939735e-07

[25]:

Draw (gfu)
Draw (gfp);

[ ]: