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.occ import *

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
Draw (shape);

geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
mesh.Curve(3)
Draw (mesh);

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

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

(u,p),(v,q) = X.TnT()

stokes = InnerProduct(Grad(u), Grad(v))*dx + div(u)*q*dx + div(v)*p*dx
a = BilinearForm(stokes).Assemble()

Set inhomogeneous Dirichlet boundary condition only on inlet boundary:

gf = GridFunction(X)
gfu, gfp = gf.components

uin = CF ( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
Draw(gfu, mesh, min=0, max=2)
SetVisualization(max=2)

Solve equation:

res = -a.mat * gf.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
gf.vec.data += inv * res
Draw(gfu, mesh);

Testing different velocity-pressure pairs

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

def SolveStokes(X):
    (u,p),(v,q) = X.TnT()

    stokes = InnerProduct(Grad(u), Grad(v))*dx + div(u)*q*dx + div(v)*p*dx
    a = BilinearForm(stokes).Assemble()

    gf = GridFunction(X)
    gfu, gfp = gf.components

    uin = CF ( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
    gfu.Set(uin, definedon=mesh.Boundaries("inlet"))

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

    Draw(gfu)

    return gfu

Higher order Taylor-Hood elements:

V = VectorH1(mesh, order=4, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=3)
X = V*Q

gfu = SolveStokes(X)

With discontinuous pressure elements P2-P1 is unstable:

V = VectorH1(mesh, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = V*Q

gfu = SolveStokes(X)

V.ndof = 3368 , Q.ndof = 2364

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Cell In[9], line 6
      3 print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
      4 X = V*Q
----> 6 gfu = SolveStokes(X)

Cell In[7], line 14, in SolveStokes(X)
     11 gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
     13 res = -a.mat * gf.vec
---> 14 inv = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack")
     15 gf.vec.data += inv * res
     17 Draw(gfu)

NgException: UmfpackInverse: Numeric factorization failed.

However, P2-P1 elements work on sub-divided simplicial meshes using Alfeld splits:

mesh2 = Mesh(geo.GenerateMesh(maxh=0.05)).Curve(3)
mesh2.SplitElements_Alfeld()
V = VectorH1(mesh2, order=2, dirichlet="wall|inlet|cyl")
Q = L2(mesh2, order=1)
print ("V.ndof =", V.ndof, ", Q.ndof =", Q.ndof)
X = V*Q
gfu = SolveStokes(X)

V.ndof = 9672 , Q.ndof = 7092

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

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

gfu = SolveStokes(X)

4944
V.ndof = 4944 , Q.ndof = 2364

the mini element:

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

gfu = SolveStokes(X)

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(InnerProduct(Grad(u),Grad(v))*dx).Assemble()
b = BilinearForm(div(u)*q*dx).Assemble()

Needed as preconditioner for the pressure:

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

Two right hand sides for the two spaces:

f = LinearForm(V).Assemble()
g = LinearForm(Q).Assemble();

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, printrates=True, initialize=False);

LinearSolver iteration 1, residual = 4.608947775183432
LinearSolver iteration 2, residual = 2.3205138844581272
LinearSolver iteration 3, residual = 1.9304918812163188
LinearSolver iteration 4, residual = 1.5339634630760308
LinearSolver iteration 5, residual = 1.486615746949375
LinearSolver iteration 6, residual = 1.2774281643545382
LinearSolver iteration 7, residual = 1.2336050111420545
LinearSolver iteration 8, residual = 1.0767167802433224
LinearSolver iteration 9, residual = 0.9764208183923561
LinearSolver iteration 10, residual = 0.8124102954934218
LinearSolver iteration 11, residual = 0.7630326288071915
LinearSolver iteration 12, residual = 0.6710519572747312
LinearSolver iteration 13, residual = 0.6468552815219277
LinearSolver iteration 14, residual = 0.5965412205029288
LinearSolver iteration 15, residual = 0.5729173657364452
LinearSolver iteration 16, residual = 0.5322237098138067
LinearSolver iteration 17, residual = 0.5244109437623008
LinearSolver iteration 18, residual = 0.49120023781930605
LinearSolver iteration 19, residual = 0.4816302555990854
LinearSolver iteration 20, residual = 0.4549592480777766
LinearSolver iteration 21, residual = 0.43538867555473565
LinearSolver iteration 22, residual = 0.3922054493950564
LinearSolver iteration 23, residual = 0.3702778832189983
LinearSolver iteration 24, residual = 0.3011003348307074
LinearSolver iteration 25, residual = 0.28920675378549193
LinearSolver iteration 26, residual = 0.17308383727234683
LinearSolver iteration 27, residual = 0.16928179934914414
LinearSolver iteration 28, residual = 0.0865744173194065
LinearSolver iteration 29, residual = 0.08562898863020547
LinearSolver iteration 30, residual = 0.04685445234258132
LinearSolver iteration 31, residual = 0.04660020674567219
LinearSolver iteration 32, residual = 0.023441789137753882
LinearSolver iteration 33, residual = 0.02011841137302586
LinearSolver iteration 34, residual = 0.013048209565297984
LinearSolver iteration 35, residual = 0.012480171506022307
LinearSolver iteration 36, residual = 0.00913530597535095
LinearSolver iteration 37, residual = 0.008779885649182962
LinearSolver iteration 38, residual = 0.006460732446737465
LinearSolver iteration 39, residual = 0.006058729455791614
LinearSolver iteration 40, residual = 0.004153069889943506
LinearSolver iteration 41, residual = 0.003977244990867536
LinearSolver iteration 42, residual = 0.0025498351992947602
LinearSolver iteration 43, residual = 0.0023629754253478344
LinearSolver iteration 44, residual = 0.0014897024969475145
LinearSolver iteration 45, residual = 0.0014380081708648003
LinearSolver iteration 46, residual = 0.0009456247234569049
LinearSolver iteration 47, residual = 0.0009042255078697641
LinearSolver iteration 48, residual = 0.0006050248304508098
LinearSolver iteration 49, residual = 0.000565021621693655
LinearSolver iteration 50, residual = 0.0003967204344210622
LinearSolver iteration 51, residual = 0.00037103618432564247
LinearSolver iteration 52, residual = 0.00024447080287516224
LinearSolver iteration 53, residual = 0.00023962944985237472
LinearSolver iteration 54, residual = 0.00013810579331519104
LinearSolver iteration 55, residual = 0.00012607029920635266
LinearSolver iteration 56, residual = 6.197735441205766e-05
LinearSolver iteration 57, residual = 5.419189162627352e-05
LinearSolver iteration 58, residual = 2.681788263677865e-05
LinearSolver iteration 59, residual = 2.3018259029950144e-05
LinearSolver iteration 60, residual = 1.0743945032630918e-05
LinearSolver iteration 61, residual = 1.0198694986631752e-05
LinearSolver iteration 62, residual = 4.665057163737822e-06
LinearSolver iteration 63, residual = 4.369820590533772e-06
LinearSolver iteration 64, residual = 1.8770130715748773e-06
LinearSolver iteration 65, residual = 1.5246977870812545e-06
LinearSolver iteration 66, residual = 6.167769748719773e-07
LinearSolver iteration 67, residual = 5.436025690541214e-07
LinearSolver iteration 68, residual = 2.5093238287018876e-07

Draw (gfu);

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(p*q*dx).Assemble()

f = LinearForm(V)
f.Assemble()

g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")
uin = CF( (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, printrates=True, initialize=False, maxsteps=200);

LinearSolver iteration 1, residual = 7.000179411319217
LinearSolver iteration 2, residual = 3.083435005446363
LinearSolver iteration 3, residual = 2.1164861496204947
LinearSolver iteration 4, residual = 1.8366080623996377
LinearSolver iteration 5, residual = 1.5895645511165308
LinearSolver iteration 6, residual = 1.5441710885922555
LinearSolver iteration 7, residual = 1.4644985020244865
LinearSolver iteration 8, residual = 1.3521655569577122
LinearSolver iteration 9, residual = 1.3361032449540493
LinearSolver iteration 10, residual = 1.2595625442857048
LinearSolver iteration 11, residual = 1.1495024170163213
LinearSolver iteration 12, residual = 1.1491221325036693
LinearSolver iteration 13, residual = 1.058358285227885
LinearSolver iteration 14, residual = 0.9773598464388686
LinearSolver iteration 15, residual = 0.9756045310923616
LinearSolver iteration 16, residual = 0.8670011646769753
LinearSolver iteration 17, residual = 0.7989406664024945
LinearSolver iteration 18, residual = 0.7843579390403254
LinearSolver iteration 19, residual = 0.6959411722935329
LinearSolver iteration 20, residual = 0.6804255085861403
LinearSolver iteration 21, residual = 0.6616381813443353
LinearSolver iteration 22, residual = 0.6248498919661596
LinearSolver iteration 23, residual = 0.6191429972376716
LinearSolver iteration 24, residual = 0.6067983529661147
LinearSolver iteration 25, residual = 0.5794521829418534
LinearSolver iteration 26, residual = 0.5762377937624361
LinearSolver iteration 27, residual = 0.5583795182268426
LinearSolver iteration 28, residual = 0.5344331575565257
LinearSolver iteration 29, residual = 0.5342817943561798
LinearSolver iteration 30, residual = 0.5139643507350379
LinearSolver iteration 31, residual = 0.4969832846104758
LinearSolver iteration 32, residual = 0.495959100166969
LinearSolver iteration 33, residual = 0.47682299150724156
LinearSolver iteration 34, residual = 0.4730163146750837
LinearSolver iteration 35, residual = 0.4616152701668645
LinearSolver iteration 36, residual = 0.4418922148215052
LinearSolver iteration 37, residual = 0.44184670691877065
LinearSolver iteration 38, residual = 0.4212205278983125
LinearSolver iteration 39, residual = 0.3993819042578297
LinearSolver iteration 40, residual = 0.39837889704076945
LinearSolver iteration 41, residual = 0.363098211653435
LinearSolver iteration 42, residual = 0.3484124801238035
LinearSolver iteration 43, residual = 0.32640008128436565
LinearSolver iteration 44, residual = 0.2764309783815207
LinearSolver iteration 45, residual = 0.2764258850664106
LinearSolver iteration 46, residual = 0.22603268049143824
LinearSolver iteration 47, residual = 0.19615345110067078
LinearSolver iteration 48, residual = 0.18637317691621924
LinearSolver iteration 49, residual = 0.13823609788156507
LinearSolver iteration 50, residual = 0.13250126874035537
LinearSolver iteration 51, residual = 0.11243751436660657
LinearSolver iteration 52, residual = 0.08207853067973471
LinearSolver iteration 53, residual = 0.08173368881244719
LinearSolver iteration 54, residual = 0.06039735763812966
LinearSolver iteration 55, residual = 0.04495603845786374
LinearSolver iteration 56, residual = 0.04333335843703279
LinearSolver iteration 57, residual = 0.026912433545662313
LinearSolver iteration 58, residual = 0.023795887364582075
LinearSolver iteration 59, residual = 0.019506048626305797
LinearSolver iteration 60, residual = 0.0133359745095983
LinearSolver iteration 61, residual = 0.013253249776554307
LinearSolver iteration 62, residual = 0.009730686991247926
LinearSolver iteration 63, residual = 0.007615282174447439
LinearSolver iteration 64, residual = 0.007201992059008289
LinearSolver iteration 65, residual = 0.005138225366806596
LinearSolver iteration 66, residual = 0.005061934778752756
LinearSolver iteration 67, residual = 0.0042383527900331465
LinearSolver iteration 68, residual = 0.0037078024600000017
LinearSolver iteration 69, residual = 0.003631483458203894
LinearSolver iteration 70, residual = 0.0028329868150136803
LinearSolver iteration 71, residual = 0.0026378867952036285
LinearSolver iteration 72, residual = 0.0021883492388023247
LinearSolver iteration 73, residual = 0.0015223538715951302
LinearSolver iteration 74, residual = 0.0015223536072792908
LinearSolver iteration 75, residual = 0.0010546009855633219
LinearSolver iteration 76, residual = 0.0008252543088433786
LinearSolver iteration 77, residual = 0.0007548709839584794
LinearSolver iteration 78, residual = 0.000463239479593169
LinearSolver iteration 79, residual = 0.0004556178744578989
LinearSolver iteration 80, residual = 0.00030264974868472803
LinearSolver iteration 81, residual = 0.0002380077833742326
LinearSolver iteration 82, residual = 0.0002188020335111653
LinearSolver iteration 83, residual = 0.00014362172751952203
LinearSolver iteration 84, residual = 0.00013117832451954014
LinearSolver iteration 85, residual = 0.0001089933768747302
LinearSolver iteration 86, residual = 6.922560261768888e-05
LinearSolver iteration 87, residual = 6.915158745281448e-05
LinearSolver iteration 88, residual = 4.489239962403809e-05
LinearSolver iteration 89, residual = 3.66745863125787e-05
LinearSolver iteration 90, residual = 3.156497541861379e-05
LinearSolver iteration 91, residual = 2.0186943902932957e-05
LinearSolver iteration 92, residual = 1.9643654504785657e-05
LinearSolver iteration 93, residual = 1.3523249354476612e-05
LinearSolver iteration 94, residual = 9.02921278865358e-06
LinearSolver iteration 95, residual = 8.616082652154183e-06
LinearSolver iteration 96, residual = 4.885128583211393e-06
LinearSolver iteration 97, residual = 4.2323140672951465e-06
LinearSolver iteration 98, residual = 3.2553481613098634e-06
LinearSolver iteration 99, residual = 1.929110466603289e-06
LinearSolver iteration 100, residual = 1.92492350573699e-06
LinearSolver iteration 101, residual = 1.1575115407525552e-06
LinearSolver iteration 102, residual = 8.292938314442221e-07
LinearSolver iteration 103, residual = 7.27924043457081e-07
LinearSolver iteration 104, residual = 3.9763569554595716e-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(p*q*dx).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.Assemble()

g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")
uin = CF( (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=True, initialize=False, maxsteps=200);

LinearSolver iteration 1, residual = 10.072431117077521
LinearSolver iteration 2, residual = 3.509846326863495
LinearSolver iteration 3, residual = 2.4713753452091187
LinearSolver iteration 4, residual = 1.600041654102234
LinearSolver iteration 5, residual = 1.5999941924014385
LinearSolver iteration 6, residual = 1.331848127071805
LinearSolver iteration 7, residual = 1.2842023741344113
LinearSolver iteration 8, residual = 1.0972672160421808
LinearSolver iteration 9, residual = 1.0536946299066983
LinearSolver iteration 10, residual = 0.8591400365732078
LinearSolver iteration 11, residual = 0.8019583070152705
LinearSolver iteration 12, residual = 0.684741960892436
LinearSolver iteration 13, residual = 0.6763718338306878
LinearSolver iteration 14, residual = 0.6146066572025399
LinearSolver iteration 15, residual = 0.6023824596733206
LinearSolver iteration 16, residual = 0.5555606514054319
LinearSolver iteration 17, residual = 0.5430937217136217
LinearSolver iteration 18, residual = 0.5094093630083454
LinearSolver iteration 19, residual = 0.49958044246587296
LinearSolver iteration 20, residual = 0.4687233539170855
LinearSolver iteration 21, residual = 0.4561019247305962
LinearSolver iteration 22, residual = 0.4245953748459189
LinearSolver iteration 23, residual = 0.40438140563842356
LinearSolver iteration 24, residual = 0.34555465092342474
LinearSolver iteration 25, residual = 0.3189018314414242
LinearSolver iteration 26, residual = 0.24053411996462457
LinearSolver iteration 27, residual = 0.23388758308235236
LinearSolver iteration 28, residual = 0.1319453589801984
LinearSolver iteration 29, residual = 0.12123099213408277
LinearSolver iteration 30, residual = 0.06382557163322702
LinearSolver iteration 31, residual = 0.06303180309279519
LinearSolver iteration 32, residual = 0.0308499456999827
LinearSolver iteration 33, residual = 0.029346981714250504
LinearSolver iteration 34, residual = 0.015627485603868757
LinearSolver iteration 35, residual = 0.014617906806415
LinearSolver iteration 36, residual = 0.008292871905731648
LinearSolver iteration 37, residual = 0.008150714306777567
LinearSolver iteration 38, residual = 0.0056479946169437035
LinearSolver iteration 39, residual = 0.00562076631964996
LinearSolver iteration 40, residual = 0.003804289007235055
LinearSolver iteration 41, residual = 0.003801474920013202
LinearSolver iteration 42, residual = 0.0020911677838965042
LinearSolver iteration 43, residual = 0.0020681554195363158
LinearSolver iteration 44, residual = 0.0010343859929530456
LinearSolver iteration 45, residual = 0.0009450837316212568
LinearSolver iteration 46, residual = 0.0005040197865835285
LinearSolver iteration 47, residual = 0.0004947091690999611
LinearSolver iteration 48, residual = 0.00026577591276723446
LinearSolver iteration 49, residual = 0.0002625558560124446
LinearSolver iteration 50, residual = 0.00015169641651782155
LinearSolver iteration 51, residual = 0.00015046497706181793
LinearSolver iteration 52, residual = 9.099101333531106e-05
LinearSolver iteration 53, residual = 9.064467690437921e-05
LinearSolver iteration 54, residual = 5.549445645768736e-05
LinearSolver iteration 55, residual = 5.4425552921804786e-05
LinearSolver iteration 56, residual = 3.3740357848835876e-05
LinearSolver iteration 57, residual = 3.316553714682566e-05
LinearSolver iteration 58, residual = 1.8747435764205758e-05
LinearSolver iteration 59, residual = 1.8611134264359517e-05
LinearSolver iteration 60, residual = 9.848616102008858e-06
LinearSolver iteration 61, residual = 9.695503799560197e-06
LinearSolver iteration 62, residual = 5.1974018726699e-06
LinearSolver iteration 63, residual = 5.183534411845862e-06
LinearSolver iteration 64, residual = 2.365090968966103e-06
LinearSolver iteration 65, residual = 2.301651524107217e-06
LinearSolver iteration 66, residual = 8.885736555407656e-07

Draw (gfu)
Draw (gfp);