# 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 *
import netgen.gui

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:

[2]:

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.

[3]:

ux,uy,p = X.TrialFunction()
vx,vy,q = X.TestFunction()

a = BilinearForm(X)
a.Assemble()

[3]:

<ngsolve.comp.BilinearForm at 0x7f6fe5bbf370>


Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

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

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

[6]:

def SolveStokes(X):
ux,uy,p = X.TrialFunction()
vx,vy,q = X.TestFunction()
a = BilinearForm(X)
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:

[7]:

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:

[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 = V*V*Q

gfu = SolveStokes(X)

V.ndof = 1660 , Q.ndof = 2328

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
<ipython-input-1-5ecc0da3557f> in <module>
4 X = V*V*Q
5
----> 6 gfu = SolveStokes(X)

<ipython-input-1-abe06cca85c3> 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

NgException: UmfpackInverse: Numeric factorization failed.


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

[9]:

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 = 2436 , Q.ndof = 2328


the mini element:

[10]:

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.

[11]:

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

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

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

a.Assemble()
b.Assemble()

[12]:

<ngsolve.comp.BilinearForm at 0x7f6ff2ad5530>


Needed as preconditioner for the pressure:

[13]:

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

[13]:

<ngsolve.comp.BilinearForm at 0x7f6fe5be2270>


Two right hand sides for the two spaces:

[14]:

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

g = LinearForm(Q)
g.Assemble()

[14]:

<ngsolve.comp.LinearForm at 0x7f6fb0049970>


Two GridFunctions for velocity and pressure:

[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.

[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.330331719284448
it =  1  err =  2.263376283463758
it =  2  err =  1.892595498495799
it =  3  err =  1.5277235488222025
it =  4  err =  1.4772437563202756
it =  5  err =  1.2747398116652078
it =  6  err =  1.2273571122716929
it =  7  err =  1.0796232252707203
it =  8  err =  0.9926664036457307
it =  9  err =  0.8697538447370408
it =  10  err =  0.8162079324261495
it =  11  err =  0.7259022997317126
it =  12  err =  0.7046131437605081
it =  13  err =  0.643631001630437
it =  14  err =  0.6266239595599093
it =  15  err =  0.591685129375076
it =  16  err =  0.5778780647646606
it =  17  err =  0.5491851877285633
it =  18  err =  0.532115011082765
it =  19  err =  0.5051337139491419
it =  20  err =  0.49384526542071455
it =  21  err =  0.4774921377058554
it =  22  err =  0.45999935874433057
it =  23  err =  0.43649884689453894
it =  24  err =  0.413379175507607
it =  25  err =  0.37941719659836165
it =  26  err =  0.3467393998530268
it =  27  err =  0.3051966386821344
it =  28  err =  0.2562525224248321
it =  29  err =  0.19300080894165858
it =  30  err =  0.1482938430502556
it =  31  err =  0.09777687507769889
it =  32  err =  0.0826686992149031
it =  33  err =  0.047911262763998035
it =  34  err =  0.046502023480814565
it =  35  err =  0.024526833323920062
it =  36  err =  0.02426719424161278
it =  37  err =  0.013783639910734346
it =  38  err =  0.013778935988601962
it =  39  err =  0.009300374509216055
it =  40  err =  0.009198660475496428
it =  41  err =  0.00703294426771334
it =  42  err =  0.006969996015614246
it =  43  err =  0.004351125368661292
it =  44  err =  0.004335298641310061
it =  45  err =  0.002611880528054009
it =  46  err =  0.0026044567370714297
it =  47  err =  0.0016436425856950214
it =  48  err =  0.0016162065941244563
it =  49  err =  0.0010771750925515542
it =  50  err =  0.0010514431988425407
it =  51  err =  0.00059962236355928
it =  52  err =  0.0005976951514385686
it =  53  err =  0.0004050240494687218
it =  54  err =  0.00040242226331075725
it =  55  err =  0.00025503750267261746
it =  56  err =  0.0002534442461167978
it =  57  err =  0.0001358191479420347
it =  58  err =  0.0001341113952354638
it =  59  err =  5.161676617964784e-05
it =  60  err =  5.142197738847716e-05
it =  61  err =  2.4354833745927838e-05
it =  62  err =  2.4124112714916188e-05
it =  63  err =  1.1179007810924205e-05
it =  64  err =  1.1174147344564869e-05
it =  65  err =  4.21869498205886e-06
it =  66  err =  4.2145437328840685e-06
it =  67  err =  1.8420124075366427e-06
it =  68  err =  1.8419951602022333e-06
it =  69  err =  6.566674225336547e-07
it =  70  err =  6.565843076533915e-07
it =  71  err =  2.0416947196036209e-07
it =  72  err =  2.039890739112055e-07
it =  73  err =  9.500196852887388e-08

[16]:

basevector

[17]:

Draw (gfu)

[ ]:



[ ]: