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

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 = 1664 , Q.ndof = 2334

---------------------------------------------------------------------------
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 = 2442 , Q.ndof = 2334

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

Needed as preconditioner for the pressure:

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

<ngsolve.comp.BilinearForm at 0x7fc4f79ea070>

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

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.567113174705788
LinearSolver iteration 2, residual = 2.3131306544242323
LinearSolver iteration 3, residual = 1.9197576117165607
LinearSolver iteration 4, residual = 1.5329264792077248
LinearSolver iteration 5, residual = 1.48543827450333
LinearSolver iteration 6, residual = 1.2783106453405675
LinearSolver iteration 7, residual = 1.2334349817866066
LinearSolver iteration 8, residual = 1.0824282246213506
LinearSolver iteration 9, residual = 0.995750924590334
LinearSolver iteration 10, residual = 0.8703143846374614
LinearSolver iteration 11, residual = 0.8171377680888388
LinearSolver iteration 12, residual = 0.7269213706513465
LinearSolver iteration 13, residual = 0.706207221729229
LinearSolver iteration 14, residual = 0.6440190256873731
LinearSolver iteration 15, residual = 0.6268764946004171
LinearSolver iteration 16, residual = 0.5917213405973824
LinearSolver iteration 17, residual = 0.5780662445654594
LinearSolver iteration 18, residual = 0.5493289072094152
LinearSolver iteration 19, residual = 0.5324039032206915
LinearSolver iteration 20, residual = 0.5052471794033446
LinearSolver iteration 21, residual = 0.4940233674385707
LinearSolver iteration 22, residual = 0.4775138626691476
LinearSolver iteration 23, residual = 0.46042831325251066
LinearSolver iteration 24, residual = 0.43669654174643346
LinearSolver iteration 25, residual = 0.4134659765859832
LinearSolver iteration 26, residual = 0.37916315943871476
LinearSolver iteration 27, residual = 0.34713390430404173
LinearSolver iteration 28, residual = 0.3052471056000025
LinearSolver iteration 29, residual = 0.2565800038676399
LinearSolver iteration 30, residual = 0.19161286393115326
LinearSolver iteration 31, residual = 0.14977050175884263
LinearSolver iteration 32, residual = 0.09807999310750135
LinearSolver iteration 33, residual = 0.08259027584460953
LinearSolver iteration 34, residual = 0.04732357043003127
LinearSolver iteration 35, residual = 0.046245141259990435
LinearSolver iteration 36, residual = 0.024604866245844917
LinearSolver iteration 37, residual = 0.024337515747360358
LinearSolver iteration 38, residual = 0.013922111781146086
LinearSolver iteration 39, residual = 0.013920784173787181
LinearSolver iteration 40, residual = 0.009562415358374389
LinearSolver iteration 41, residual = 0.009428591408873454
LinearSolver iteration 42, residual = 0.007150714585376771
LinearSolver iteration 43, residual = 0.007098687832640237
LinearSolver iteration 44, residual = 0.004446629197811791
LinearSolver iteration 45, residual = 0.004412245271741807
LinearSolver iteration 46, residual = 0.0025975433837623738
LinearSolver iteration 47, residual = 0.0025921631587572853
LinearSolver iteration 48, residual = 0.0016894640885924124
LinearSolver iteration 49, residual = 0.001649302458554585
LinearSolver iteration 50, residual = 0.0010653711198279352
LinearSolver iteration 51, residual = 0.001049103113046462
LinearSolver iteration 52, residual = 0.0006174817591898052
LinearSolver iteration 53, residual = 0.0006147544635181265
LinearSolver iteration 54, residual = 0.00041375056119163017
LinearSolver iteration 55, residual = 0.0004117487931243662
LinearSolver iteration 56, residual = 0.0002578161579849879
LinearSolver iteration 57, residual = 0.00025509036919492374
LinearSolver iteration 58, residual = 0.00013238071846776123
LinearSolver iteration 59, residual = 0.00013109824861237197
LinearSolver iteration 60, residual = 5.231437189439247e-05
LinearSolver iteration 61, residual = 5.2010590849767385e-05
LinearSolver iteration 62, residual = 2.4925287017848086e-05
LinearSolver iteration 63, residual = 2.4657659348114557e-05
LinearSolver iteration 64, residual = 1.0905653096273308e-05
LinearSolver iteration 65, residual = 1.0875926872751356e-05
LinearSolver iteration 66, residual = 4.194117622782712e-06
LinearSolver iteration 67, residual = 4.185367608267532e-06
LinearSolver iteration 68, residual = 1.824460679274315e-06
LinearSolver iteration 69, residual = 1.824449806889949e-06
LinearSolver iteration 70, residual = 6.228881113741227e-07
LinearSolver iteration 71, residual = 6.227781818155741e-07
LinearSolver iteration 72, residual = 1.9482176077318115e-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.946912445874184e-07

Draw (gfu)
Draw (gfp);