This page was generated from jupyter-files/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:

In [1]:

from ngsolve import *
import netgen.gui
%gui tk

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)

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

In [2]:

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

In [3]:

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 += SymbolicBFI(grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p)
a.Assemble()

Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

In [4]:

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")
Draw(Norm(velocity), mesh, "|vel|")
SetVisualization(max=2)

Solve equation:

In [5]:

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

Testing different velocity-pressure pairs¶

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

In [6]:

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 += SymbolicBFI(grad(ux)*grad(vx)+grad(uy)*grad(vy) + div_u*q + div_v*p)
    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:

In [7]:

V = H1(mesh, order=4, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=3)
X = FESpace([V,V,Q])

gfu = SolveStokes(X)

With discontinuous pressure elements P2-P1 is unstable:

In [8]:

V = H1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = FESpace([V,V,Q])

gfu = SolveStokes(X)

V.ndof = 1702 , Q.ndof = 2382

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-8-4e71f02dc134> in <module>()
      4 X = FESpace([V,V,Q])
      5
----> 6 gfu = SolveStokes(X)

<ipython-input-6-13793926ffa5> 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
     16

RuntimeError: UmfpackInverse: Numeric factorization failed.

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

In [9]:

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

gfu = SolveStokes(X)

V.ndof = 2496 , Q.ndof = 2382

the mini element:

In [10]:

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

gfu = SolveStokes(X)

VectorH1¶

A vector-valued \(H^1\)-space: Less to type and more possibilities to explore structure and optimize.

In [11]:

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

u,p = X.TrialFunction()
v,q = X.TestFunction()

a = BilinearForm(X)
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p)
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:

In [12]:

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 += SymbolicBFI(InnerProduct(grad(u),grad(v)))

b = BilinearForm(trialspace=V, testspace=Q)
b += SymbolicBFI(div(u)*q)

a.Assemble()
b.Assemble()

Needed as preconditioner for the pressure:

In [13]:

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

Two right hand sides for the two spaces:

In [14]:

f = LinearForm(V)
f += SymbolicLFI( CoefficientFunction( (0,x-0.5)) * v)
f.Assemble()

g = LinearForm(Q)
g.Assemble()

Two GridFunctions for velocity and pressure:

In [15]:

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.

In [16]:

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)

it =  0  err =  4.577250004746686
it =  1  err =  2.3137027107825414
it =  2  err =  1.9474651703872594
it =  3  err =  1.5362811821422817
it =  4  err =  1.490898743198753
it =  5  err =  1.277851941805208
it =  6  err =  1.236894318332737
it =  7  err =  1.0832164597336917
it =  8  err =  0.9954867877973922
it =  9  err =  0.8694764391752701
it =  10  err =  0.8164294710037596
it =  11  err =  0.7270898498126048
it =  12  err =  0.7062300357748325
it =  13  err =  0.6435150710145074
it =  14  err =  0.6260336874381139
it =  15  err =  0.591910771493464
it =  16  err =  0.577362791218001
it =  17  err =  0.5490895491963667
it =  18  err =  0.5299584240360132
it =  19  err =  0.504587770754635
it =  20  err =  0.4945747930809748
it =  21  err =  0.478452427389186
it =  22  err =  0.46049192930432753
it =  23  err =  0.4367980978458769
it =  24  err =  0.41255367261726367
it =  25  err =  0.38040899835370073
it =  26  err =  0.3449005478686267
it =  27  err =  0.2985421196154886
it =  28  err =  0.2527173796456808
it =  29  err =  0.17951963386513986
it =  30  err =  0.14701867167188856
it =  31  err =  0.10086316225427124
it =  32  err =  0.08594125971647408
it =  33  err =  0.049535710323364915
it =  34  err =  0.04880618304699326
it =  35  err =  0.025962616020410965
it =  36  err =  0.025723086935984515
it =  37  err =  0.01628165309575614
it =  38  err =  0.0162502272287645
it =  39  err =  0.011879675741016503
it =  40  err =  0.01145912264948227
it =  41  err =  0.0077088148705381325
it =  42  err =  0.007669232876181972
it =  43  err =  0.0047624380822534566
it =  44  err =  0.0047402737522261665
it =  45  err =  0.0024617452306963774
it =  46  err =  0.002444954453502464
it =  47  err =  0.0015421225878978428
it =  48  err =  0.0014860456371280366
it =  49  err =  0.0010744692033775487
it =  50  err =  0.0010244043188852743
it =  51  err =  0.0006385557996683751
it =  52  err =  0.0006291811667503349
it =  53  err =  0.0002990274926769061
it =  54  err =  0.0002988061757828301
it =  55  err =  0.00011824139621187269
it =  56  err =  0.00011462936306193394
it =  57  err =  5.755195930067076e-05
it =  58  err =  5.56951683872658e-05
it =  59  err =  2.4705121755237746e-05
it =  60  err =  2.470480085269251e-05
it =  61  err =  1.6824109227260327e-05
it =  62  err =  1.6572422543775576e-05
it =  63  err =  1.316326692504858e-05
it =  64  err =  1.2897753939698155e-05
it =  65  err =  9.293568917613104e-06
it =  66  err =  9.15575741204158e-06
it =  67  err =  4.022673678037305e-06
it =  68  err =  4.021710787407399e-06
it =  69  err =  1.3784244536605242e-06
it =  70  err =  1.377553514553071e-06
it =  71  err =  6.003552147432753e-07
it =  72  err =  5.882429968010617e-07
it =  73  err =  2.4024447119496915e-07
it =  74  err =  2.3576449515661013e-07
it =  75  err =  1.0050977490275618e-07
it =  76  err =  9.69663021570943e-08

Out[16]:

basevector

In [17]:

Draw (gfu)

In [ ]: