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)

 Generate Mesh from spline geometry
 Boundary mesh done, np = 108
 CalcLocalH: 108 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters
 Curve elements, order = 3
 Update clusters
 Curving edges
 Curving faces

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

assemble VOL element 774/774

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)

setvalues element 108/108

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

call umfpack ... done

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)

assemble VOL element 774/774
setvalues element 108/108
call umfpack ... done

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 = 1656 , Q.ndof = 2322
assemble VOL element 774/774
setvalues element 108/108
call umfpack ...

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
/tmp/ipykernel_28030/1189186610.py in <module>
      4 X = V*V*Q
      5
----> 6 gfu = SolveStokes(X)

/tmp/ipykernel_28030/811579872.py 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:

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 = 2430 , Q.ndof = 2322

UMFPACK V5.7.4 (Feb 1, 2016): WARNING: matrix is singular

assemble VOL element 774/774
setvalues element 108/108
call umfpack ... done

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)

assemble VOL element 774/774
setvalues element 108/108
call umfpack ... done

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)

assemble VOL element 774/774
setvalues element 108/108
call umfpack ... done

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

assemble VOL element 774/774
assemble mixed bilinearform

<ngsolve.comp.BilinearForm at 0x7fc64d247970>

Needed as preconditioner for the pressure:

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

assemble VOL element 774/774

<ngsolve.comp.BilinearForm at 0x7fc67830fb70>

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

assemble VOL element 774/774

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

setvalues element 108/108

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)

call pardiso ... done
call pardiso ... done
LinearSolver iteration 1, residual = 4.337069360146072
LinearSolver iteration 2, residual = 2.264713596298657
LinearSolver iteration 3, residual = 1.892551008538201
LinearSolver iteration 4, residual = 1.5278298950439118
LinearSolver iteration 5, residual = 1.4774123933006966
LinearSolver iteration 6, residual = 1.2748646703687811
LinearSolver iteration 7, residual = 1.227350244839074
LinearSolver iteration 8, residual = 1.0796748148348154
LinearSolver iteration 9, residual = 0.9927213467063855
LinearSolver iteration 10, residual = 0.8697728301343073
LinearSolver iteration 11, residual = 0.816271182836896
LinearSolver iteration 12, residual = 0.7259673235404897
LinearSolver iteration 13, residual = 0.7046418006753331
LinearSolver iteration 14, residual = 0.6436420998832291
LinearSolver iteration 15, residual = 0.6266519183176467
LinearSolver iteration 16, residual = 0.5916941058249604
LinearSolver iteration 17, residual = 0.5778872637584644
LinearSolver iteration 18, residual = 0.5491983683896285
LinearSolver iteration 19, residual = 0.5321339542232276
LinearSolver iteration 20, residual = 0.5051520176131404
LinearSolver iteration 21, residual = 0.49388035273938846
LinearSolver iteration 22, residual = 0.4774977661080504
LinearSolver iteration 23, residual = 0.46005121716109815
LinearSolver iteration 24, residual = 0.43657068469684934
LinearSolver iteration 25, residual = 0.41356565003754386
LinearSolver iteration 26, residual = 0.3795032471424406
LinearSolver iteration 27, residual = 0.3467761879111833
LinearSolver iteration 28, residual = 0.30529859256434966
LinearSolver iteration 29, residual = 0.25676736950680873
LinearSolver iteration 30, residual = 0.1935837308698832
LinearSolver iteration 31, residual = 0.14865328514139542
LinearSolver iteration 32, residual = 0.09766275665495895
LinearSolver iteration 33, residual = 0.0825927053478491
LinearSolver iteration 34, residual = 0.04796176411133548
LinearSolver iteration 35, residual = 0.046561374422869194
LinearSolver iteration 36, residual = 0.024448642105259134
LinearSolver iteration 37, residual = 0.024187420660011078
LinearSolver iteration 38, residual = 0.013787348456802321
LinearSolver iteration 39, residual = 0.013782015309911428
LinearSolver iteration 40, residual = 0.009231401649420294
LinearSolver iteration 41, residual = 0.009136165052020364
LinearSolver iteration 42, residual = 0.006992878478949685
LinearSolver iteration 43, residual = 0.006933444411962679
LinearSolver iteration 44, residual = 0.004323990920530351
LinearSolver iteration 45, residual = 0.004307641889957301
LinearSolver iteration 46, residual = 0.002585414701387047
LinearSolver iteration 47, residual = 0.002579809247615046
LinearSolver iteration 48, residual = 0.0016354435576020616
LinearSolver iteration 49, residual = 0.0016073271640315529
LinearSolver iteration 50, residual = 0.0010691425353595155
LinearSolver iteration 51, residual = 0.0010456470058024623
LinearSolver iteration 52, residual = 0.000595759153615959
LinearSolver iteration 53, residual = 0.0005939215989897622
LinearSolver iteration 54, residual = 0.00040290386616996244
LinearSolver iteration 55, residual = 0.00040066380782038003
LinearSolver iteration 56, residual = 0.0002543005761819086
LinearSolver iteration 57, residual = 0.0002527535623158061
LinearSolver iteration 58, residual = 0.00013526924429052613
LinearSolver iteration 59, residual = 0.00013376572240051968
LinearSolver iteration 60, residual = 5.1163604302161335e-05
LinearSolver iteration 61, residual = 5.0942439178644265e-05
LinearSolver iteration 62, residual = 2.443219452273615e-05
LinearSolver iteration 63, residual = 2.4186932007682973e-05
LinearSolver iteration 64, residual = 1.113132543105798e-05
LinearSolver iteration 65, residual = 1.1128361647134322e-05
LinearSolver iteration 66, residual = 4.244148442199896e-06
LinearSolver iteration 67, residual = 4.2413246789289464e-06
LinearSolver iteration 68, residual = 1.8252703488600543e-06
LinearSolver iteration 69, residual = 1.8248987304295922e-06
LinearSolver iteration 70, residual = 6.602943285291293e-07
LinearSolver iteration 71, residual = 6.602098965276247e-07
LinearSolver iteration 72, residual = 2.049931114192133e-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()

 Generate Mesh from spline geometry
 Boundary mesh done, np = 258
 CalcLocalH: 258 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters
assemble VOL element 4548/4548
assemble mixed bilinearform
assemble VOL element 4548/4548

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

assemble VOL element 4548/4548
assemble VOL element 4548/4548
setvalues element 258/258
call pardiso ... done
call pardiso ... done
LinearSolver iteration 1, residual = 7.093524574559076
LinearSolver iteration 2, residual = 3.123301846060359
LinearSolver iteration 3, residual = 2.1396216166041753
LinearSolver iteration 4, residual = 1.914219462615039
LinearSolver iteration 5, residual = 1.6116196031976782
LinearSolver iteration 6, residual = 1.5489603152856475
LinearSolver iteration 7, residual = 1.4848591185063968
LinearSolver iteration 8, residual = 1.3800580493063581
LinearSolver iteration 9, residual = 1.3386235485379536
LinearSolver iteration 10, residual = 1.284152849911857
LinearSolver iteration 11, residual = 1.171400089134812
LinearSolver iteration 12, residual = 1.1640324250820986
LinearSolver iteration 13, residual = 1.0888266346366466
LinearSolver iteration 14, residual = 1.0126085439860515
LinearSolver iteration 15, residual = 1.0008115366910437
LinearSolver iteration 16, residual = 0.8910955738823616
LinearSolver iteration 17, residual = 0.8729880985957466
LinearSolver iteration 18, residual = 0.7981674440668808
LinearSolver iteration 19, residual = 0.7244171610865654
LinearSolver iteration 20, residual = 0.7233058879256861
LinearSolver iteration 21, residual = 0.6687530797250423
LinearSolver iteration 22, residual = 0.6526056449597718
LinearSolver iteration 23, residual = 0.6436431389122804
LinearSolver iteration 24, residual = 0.6159761414220264
LinearSolver iteration 25, residual = 0.6115222527483778
LinearSolver iteration 26, residual = 0.5991667862423437
LinearSolver iteration 27, residual = 0.5736793196432238
LinearSolver iteration 28, residual = 0.5714962016331891
LinearSolver iteration 29, residual = 0.5542631844122419
LinearSolver iteration 30, residual = 0.532055024445272
LinearSolver iteration 31, residual = 0.531921457390652
LinearSolver iteration 32, residual = 0.510729746762853
LinearSolver iteration 33, residual = 0.5019586208536521
LinearSolver iteration 34, residual = 0.49560770088990036
LinearSolver iteration 35, residual = 0.4750072960858663
LinearSolver iteration 36, residual = 0.4713090543262219
LinearSolver iteration 37, residual = 0.46146376709121284
LinearSolver iteration 38, residual = 0.4391832079345872
LinearSolver iteration 39, residual = 0.437368683586536
LinearSolver iteration 40, residual = 0.4211512919810195
LinearSolver iteration 41, residual = 0.39367115138774345
LinearSolver iteration 42, residual = 0.39347657129658536
LinearSolver iteration 43, residual = 0.36001180083278406
LinearSolver iteration 44, residual = 0.33039855672199453
LinearSolver iteration 45, residual = 0.3246152628977096
LinearSolver iteration 46, residual = 0.26709807045403816
LinearSolver iteration 47, residual = 0.25593084790674575
LinearSolver iteration 48, residual = 0.22792999737229336
LinearSolver iteration 49, residual = 0.18718749769798987
LinearSolver iteration 50, residual = 0.18567357579428906
LinearSolver iteration 51, residual = 0.1351835183435528
LinearSolver iteration 52, residual = 0.12299160854248277
LinearSolver iteration 53, residual = 0.10499773205792497
LinearSolver iteration 54, residual = 0.07415268727177486
LinearSolver iteration 55, residual = 0.07292909225916672
LinearSolver iteration 56, residual = 0.0534475562309993
LinearSolver iteration 57, residual = 0.03862119936371087
LinearSolver iteration 58, residual = 0.03670906203305091
LinearSolver iteration 59, residual = 0.023268422282318204
LinearSolver iteration 60, residual = 0.022768409341301463
LinearSolver iteration 61, residual = 0.0166865078791172
LinearSolver iteration 62, residual = 0.012517856611925551
LinearSolver iteration 63, residual = 0.012287318123687411
LinearSolver iteration 64, residual = 0.00835034915350182
LinearSolver iteration 65, residual = 0.007108799169281131
LinearSolver iteration 66, residual = 0.006255289192429071
LinearSolver iteration 67, residual = 0.004658268592092344
LinearSolver iteration 68, residual = 0.004614851095049302
LinearSolver iteration 69, residual = 0.003926067368038764
LinearSolver iteration 70, residual = 0.0033784110116021836
LinearSolver iteration 71, residual = 0.003339413114352952
LinearSolver iteration 72, residual = 0.0025256152680784155
LinearSolver iteration 73, residual = 0.0022577605196393885
LinearSolver iteration 74, residual = 0.0019389233894438346
LinearSolver iteration 75, residual = 0.0013265439371479171
LinearSolver iteration 76, residual = 0.0013047760787624745
LinearSolver iteration 77, residual = 0.0009821054137240845
LinearSolver iteration 78, residual = 0.0006879917754877875
LinearSolver iteration 79, residual = 0.0006866656904256515
LinearSolver iteration 80, residual = 0.0004627504518609887
LinearSolver iteration 81, residual = 0.0003366306396003939
LinearSolver iteration 82, residual = 0.0003247538495450309
LinearSolver iteration 83, residual = 0.00020472316694769625
LinearSolver iteration 84, residual = 0.00017302449351049687
LinearSolver iteration 85, residual = 0.00015055955198062793
LinearSolver iteration 86, residual = 9.503979483358908e-05
LinearSolver iteration 87, residual = 9.139967370081062e-05
LinearSolver iteration 88, residual = 6.205580470675907e-05
LinearSolver iteration 89, residual = 4.417482518300715e-05
LinearSolver iteration 90, residual = 4.208015931432418e-05
LinearSolver iteration 91, residual = 2.639415198204693e-05
LinearSolver iteration 92, residual = 2.3198211487033754e-05
LinearSolver iteration 93, residual = 2.05976446848021e-05
LinearSolver iteration 94, residual = 1.335750488850388e-05
LinearSolver iteration 95, residual = 1.178183346676212e-05
LinearSolver iteration 96, residual = 9.221492579476037e-06
LinearSolver iteration 97, residual = 5.430999939758362e-06
LinearSolver iteration 98, residual = 5.424479296528133e-06
LinearSolver iteration 99, residual = 3.2829817157318746e-06
LinearSolver iteration 100, residual = 2.3963406655547116e-06
LinearSolver iteration 101, residual = 2.1384549823450807e-06
LinearSolver iteration 102, residual = 1.1841165235066038e-06
LinearSolver iteration 103, residual = 1.1037166205473494e-06
LinearSolver iteration 104, residual = 7.415661381672827e-07
LinearSolver iteration 105, residual = 4.5723926977354625e-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()

 Generate Mesh from spline geometry
 Boundary mesh done, np = 258
 CalcLocalH: 258 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters
assemble VOL element 4548/4548
 Refine mesh
 Update mesh topology
 Update clusters

# 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);

assemble VOL element 18192/18192
assemble mixed bilinearform
assemble VOL element 18192/18192
setvalues element 516/516
call pardiso ... done
LinearSolver converged in 70 iterations to residual 6.656296451166658e-07

Draw (gfu)
Draw (gfp);