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)