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

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 = 1628 , Q.ndof = 2280

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Input In [8], in <cell line: 6>()
      3 print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
      4 X = V*V*Q
----> 6 gfu = SolveStokes(X)

Input In [6], 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
     18 velocity = CoefficientFunction(gfu.components[0:2])

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 = 2388 , Q.ndof = 2280

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

Needed as preconditioner for the pressure:

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

<ngsolve.comp.BilinearForm at 0x7f555f4cc7b0>

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

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.526596218952118
LinearSolver iteration 2, residual = 2.304142025958604
LinearSolver iteration 3, residual = 1.9196740519582585
LinearSolver iteration 4, residual = 1.5323644425481462
LinearSolver iteration 5, residual = 1.484826347682785
LinearSolver iteration 6, residual = 1.277601132450486
LinearSolver iteration 7, residual = 1.2329354080834554
LinearSolver iteration 8, residual = 1.0820254610603615
LinearSolver iteration 9, residual = 0.9948403735097223
LinearSolver iteration 10, residual = 0.8699950467604349
LinearSolver iteration 11, residual = 0.8170982764599389
LinearSolver iteration 12, residual = 0.7266184452514605
LinearSolver iteration 13, residual = 0.7061370158841307
LinearSolver iteration 14, residual = 0.6439587414074485
LinearSolver iteration 15, residual = 0.6266376032903321
LinearSolver iteration 16, residual = 0.5916319320633947
LinearSolver iteration 17, residual = 0.5781591071676442
LinearSolver iteration 18, residual = 0.5493006210789705
LinearSolver iteration 19, residual = 0.5322456336368061
LinearSolver iteration 20, residual = 0.5051647310005571
LinearSolver iteration 21, residual = 0.49413174083465417
LinearSolver iteration 22, residual = 0.47752994974294233
LinearSolver iteration 23, residual = 0.46061888939187023
LinearSolver iteration 24, residual = 0.4367177947588869
LinearSolver iteration 25, residual = 0.4138130447682085
LinearSolver iteration 26, residual = 0.37935164234397434
LinearSolver iteration 27, residual = 0.3477107529720553
LinearSolver iteration 28, residual = 0.305484696945418
LinearSolver iteration 29, residual = 0.256942468275652
LinearSolver iteration 30, residual = 0.19217073708063925
LinearSolver iteration 31, residual = 0.1505681067527332
LinearSolver iteration 32, residual = 0.09801077954924
LinearSolver iteration 33, residual = 0.08268533255919523
LinearSolver iteration 34, residual = 0.04736873878497305
LinearSolver iteration 35, residual = 0.04627361839693014
LinearSolver iteration 36, residual = 0.02436433598776867
LinearSolver iteration 37, residual = 0.024125148875699616
LinearSolver iteration 38, residual = 0.013670010485413114
LinearSolver iteration 39, residual = 0.013669645643420733
LinearSolver iteration 40, residual = 0.00926012926623174
LinearSolver iteration 41, residual = 0.009166550359184589
LinearSolver iteration 42, residual = 0.006917801972287729
LinearSolver iteration 43, residual = 0.0068806702936758425
LinearSolver iteration 44, residual = 0.00430049067721782
LinearSolver iteration 45, residual = 0.00427829183203173
LinearSolver iteration 46, residual = 0.002519096243535504
LinearSolver iteration 47, residual = 0.0025170536759848185
LinearSolver iteration 48, residual = 0.0016345661906415473
LinearSolver iteration 49, residual = 0.001601608348364574
LinearSolver iteration 50, residual = 0.0010386796634071667
LinearSolver iteration 51, residual = 0.0010248992269996015
LinearSolver iteration 52, residual = 0.0005974607197787437
LinearSolver iteration 53, residual = 0.0005950321898701294
LinearSolver iteration 54, residual = 0.00040453686909553864
LinearSolver iteration 55, residual = 0.00040334134723656387
LinearSolver iteration 56, residual = 0.0002516257718836788
LinearSolver iteration 57, residual = 0.00024971829558215257
LinearSolver iteration 58, residual = 0.00013038344482102786
LinearSolver iteration 59, residual = 0.00012980847957816524
LinearSolver iteration 60, residual = 5.1355352978953214e-05
LinearSolver iteration 61, residual = 5.0997881661856817e-05
LinearSolver iteration 62, residual = 2.4647726880399786e-05
LinearSolver iteration 63, residual = 2.4491353739239848e-05
LinearSolver iteration 64, residual = 1.0731546091995218e-05
LinearSolver iteration 65, residual = 1.0721374391683005e-05
LinearSolver iteration 66, residual = 4.245914588356047e-06
LinearSolver iteration 67, residual = 4.245116290118095e-06
LinearSolver iteration 68, residual = 1.8736774304270876e-06
LinearSolver iteration 69, residual = 1.871176846150075e-06
LinearSolver iteration 70, residual = 6.34118210119666e-07
LinearSolver iteration 71, residual = 6.340980264968924e-07
LinearSolver iteration 72, residual = 2.0742648325677248e-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.0847849220105745
LinearSolver iteration 2, residual = 3.1187259891329164
LinearSolver iteration 3, residual = 2.1379046608229504
LinearSolver iteration 4, residual = 1.9074206900544455
LinearSolver iteration 5, residual = 1.6100604110336454
LinearSolver iteration 6, residual = 1.5486863950378131
LinearSolver iteration 7, residual = 1.4823210375932343
LinearSolver iteration 8, residual = 1.3789544944161232
LinearSolver iteration 9, residual = 1.3384923437864902
LinearSolver iteration 10, residual = 1.2791276590573433
LinearSolver iteration 11, residual = 1.1705427675317013
LinearSolver iteration 12, residual = 1.1633499593931433
LinearSolver iteration 13, residual = 1.0859851203190911
LinearSolver iteration 14, residual = 1.0113563228852067
LinearSolver iteration 15, residual = 0.9975122564877411
LinearSolver iteration 16, residual = 0.8874929184937893
LinearSolver iteration 17, residual = 0.8715250435688432
LinearSolver iteration 18, residual = 0.7958684446339616
LinearSolver iteration 19, residual = 0.724416417153301
LinearSolver iteration 20, residual = 0.7237506544094237
LinearSolver iteration 21, residual = 0.6708548712575405
LinearSolver iteration 22, residual = 0.6522193987033609
LinearSolver iteration 23, residual = 0.6436823476113409
LinearSolver iteration 24, residual = 0.6156999688073403
LinearSolver iteration 25, residual = 0.6112238937739011
LinearSolver iteration 26, residual = 0.5987753990356403
LinearSolver iteration 27, residual = 0.5740182315186663
LinearSolver iteration 28, residual = 0.5712995187952564
LinearSolver iteration 29, residual = 0.5551596781670438
LinearSolver iteration 30, residual = 0.5326568028197286
LinearSolver iteration 31, residual = 0.5326560798952301
LinearSolver iteration 32, residual = 0.5115997017048597
LinearSolver iteration 33, residual = 0.5017459046659115
LinearSolver iteration 34, residual = 0.4969346171683325
LinearSolver iteration 35, residual = 0.4763346942718254
LinearSolver iteration 36, residual = 0.4712856839812733
LinearSolver iteration 37, residual = 0.4607108422489289
LinearSolver iteration 38, residual = 0.4383192107985963
LinearSolver iteration 39, residual = 0.4368673746736316
LinearSolver iteration 40, residual = 0.42039726478419914
LinearSolver iteration 41, residual = 0.3937475832103462
LinearSolver iteration 42, residual = 0.3933763574745416
LinearSolver iteration 43, residual = 0.36028377839096865
LinearSolver iteration 44, residual = 0.32961763056079274
LinearSolver iteration 45, residual = 0.3242421789306683
LinearSolver iteration 46, residual = 0.267297694844891
LinearSolver iteration 47, residual = 0.2550075283807847
LinearSolver iteration 48, residual = 0.22770481990030852
LinearSolver iteration 49, residual = 0.18588341149074766
LinearSolver iteration 50, residual = 0.18414883175709623
LinearSolver iteration 51, residual = 0.13314679787448797
LinearSolver iteration 52, residual = 0.12229260761704674
LinearSolver iteration 53, residual = 0.10455316317548968
LinearSolver iteration 54, residual = 0.07366698305972842
LinearSolver iteration 55, residual = 0.07251347614707696
LinearSolver iteration 56, residual = 0.05287862113540375
LinearSolver iteration 57, residual = 0.038350865367667124
LinearSolver iteration 58, residual = 0.03627053273978728
LinearSolver iteration 59, residual = 0.022972988173685652
LinearSolver iteration 60, residual = 0.022578620106767856
LinearSolver iteration 61, residual = 0.016305069975769918
LinearSolver iteration 62, residual = 0.012374107528290237
LinearSolver iteration 63, residual = 0.01203398210131206
LinearSolver iteration 64, residual = 0.00822565033775119
LinearSolver iteration 65, residual = 0.007053974275159431
LinearSolver iteration 66, residual = 0.0061886437683723065
LinearSolver iteration 67, residual = 0.004641883148217402
LinearSolver iteration 68, residual = 0.00460274936940639
LinearSolver iteration 69, residual = 0.003914378946518575
LinearSolver iteration 70, residual = 0.003380542152779559
LinearSolver iteration 71, residual = 0.0033501470853653883
LinearSolver iteration 72, residual = 0.0025584609207297687
LinearSolver iteration 73, residual = 0.002244827637426911
LinearSolver iteration 74, residual = 0.0019933321147920804
LinearSolver iteration 75, residual = 0.0013546533131414836
LinearSolver iteration 76, residual = 0.0013041633600495825
LinearSolver iteration 77, residual = 0.0009950947446727685
LinearSolver iteration 78, residual = 0.0006778518480657907
LinearSolver iteration 79, residual = 0.0006768927406417723
LinearSolver iteration 80, residual = 0.0004186126951000008
LinearSolver iteration 81, residual = 0.0003269047474549921
LinearSolver iteration 82, residual = 0.0003060562015474792
LinearSolver iteration 83, residual = 0.00020001283993460031
LinearSolver iteration 84, residual = 0.0001672569028277847
LinearSolver iteration 85, residual = 0.00014955148777428612
LinearSolver iteration 86, residual = 9.248112978694838e-05
LinearSolver iteration 87, residual = 8.851434405849625e-05
LinearSolver iteration 88, residual = 6.202971738491974e-05
LinearSolver iteration 89, residual = 4.281889604756001e-05
LinearSolver iteration 90, residual = 4.141174648167259e-05
LinearSolver iteration 91, residual = 2.5957409027657297e-05
LinearSolver iteration 92, residual = 2.2470022963720288e-05
LinearSolver iteration 93, residual = 1.8917375041877856e-05
LinearSolver iteration 94, residual = 1.1649292126144114e-05
LinearSolver iteration 95, residual = 1.136538808998595e-05
LinearSolver iteration 96, residual = 7.590183158296661e-06
LinearSolver iteration 97, residual = 5.071617884861461e-06
LinearSolver iteration 98, residual = 4.90280141221119e-06
LinearSolver iteration 99, residual = 3.031285799654416e-06
LinearSolver iteration 100, residual = 2.375335418272855e-06
LinearSolver iteration 101, residual = 2.208362065702669e-06
LinearSolver iteration 102, residual = 1.250500995681854e-06
LinearSolver iteration 103, residual = 1.089982154107014e-06
LinearSolver iteration 104, residual = 8.116108773739261e-07
LinearSolver iteration 105, residual = 4.609232927264937e-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 5.946912445878173e-07

Draw (gfu)
Draw (gfp);