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.9304918812163194
LinearSolver iteration 4, residual = 1.533963463076032
LinearSolver iteration 5, residual = 1.4866157469493764
LinearSolver iteration 6, residual = 1.27742816435454
LinearSolver iteration 7, residual = 1.2336050111420565
LinearSolver iteration 8, residual = 1.076716780243325
LinearSolver iteration 9, residual = 0.9764208183923596
LinearSolver iteration 10, residual = 0.8124102954934254
LinearSolver iteration 11, residual = 0.7630326288071948
LinearSolver iteration 12, residual = 0.6710519572747341
LinearSolver iteration 13, residual = 0.6468552815219305
LinearSolver iteration 14, residual = 0.5965412205029312
LinearSolver iteration 15, residual = 0.5729173657364472
LinearSolver iteration 16, residual = 0.5322237098138083
LinearSolver iteration 17, residual = 0.5244109437623018
LinearSolver iteration 18, residual = 0.49120023781930633
LinearSolver iteration 19, residual = 0.48163025559908507
LinearSolver iteration 20, residual = 0.45495924807777527
LinearSolver iteration 21, residual = 0.4353886755547329
LinearSolver iteration 22, residual = 0.3922054493950513
LinearSolver iteration 23, residual = 0.3702778832189908
LinearSolver iteration 24, residual = 0.3011003348306984
LinearSolver iteration 25, residual = 0.28920675378548133
LinearSolver iteration 26, residual = 0.17308383727233664
LinearSolver iteration 27, residual = 0.16928179934913418
LinearSolver iteration 28, residual = 0.08657441731939673
LinearSolver iteration 29, residual = 0.08562898863019655
LinearSolver iteration 30, residual = 0.046854452342577256
LinearSolver iteration 31, residual = 0.04660020674566824
LinearSolver iteration 32, residual = 0.023441789137752054
LinearSolver iteration 33, residual = 0.020118411373024562
LinearSolver iteration 34, residual = 0.013048209565296945
LinearSolver iteration 35, residual = 0.012480171506021435
LinearSolver iteration 36, residual = 0.009135305975350208
LinearSolver iteration 37, residual = 0.008779885649182239
LinearSolver iteration 38, residual = 0.006460732446736903
LinearSolver iteration 39, residual = 0.006058729455791058
LinearSolver iteration 40, residual = 0.0041530698899431276
LinearSolver iteration 41, residual = 0.0039772449908671735
LinearSolver iteration 42, residual = 0.0025498351992945057
LinearSolver iteration 43, residual = 0.0023629754253476232
LinearSolver iteration 44, residual = 0.001489702496947371
LinearSolver iteration 45, residual = 0.001438008170864664
LinearSolver iteration 46, residual = 0.0009456247234568032
LinearSolver iteration 47, residual = 0.0009042255078696641
LinearSolver iteration 48, residual = 0.0006050248304507476
LinearSolver iteration 49, residual = 0.0005650216216935938
LinearSolver iteration 50, residual = 0.0003967204344210219
LinearSolver iteration 51, residual = 0.00037103618432560674
LinearSolver iteration 52, residual = 0.00024447080287513654
LinearSolver iteration 53, residual = 0.00023962944985235049
LinearSolver iteration 54, residual = 0.00013810579331517543
LinearSolver iteration 55, residual = 0.0001260702992063386
LinearSolver iteration 56, residual = 6.197735441204934e-05
LinearSolver iteration 57, residual = 5.4191891626266326e-05
LinearSolver iteration 58, residual = 2.681788263677515e-05
LinearSolver iteration 59, residual = 2.3018259029946735e-05
LinearSolver iteration 60, residual = 1.0743945032629388e-05
LinearSolver iteration 61, residual = 1.0198694986630182e-05
LinearSolver iteration 62, residual = 4.665057163737225e-06
LinearSolver iteration 63, residual = 4.369820590533216e-06
LinearSolver iteration 64, residual = 1.8770130715746128e-06
LinearSolver iteration 65, residual = 1.5246977870810158e-06
LinearSolver iteration 66, residual = 6.167769748718805e-07
LinearSolver iteration 67, residual = 5.436025690540095e-07
LinearSolver iteration 68, residual = 2.5093238287015657e-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.116486149620493
LinearSolver iteration 4, residual = 1.8366080623996368
LinearSolver iteration 5, residual = 1.5895645511165304
LinearSolver iteration 6, residual = 1.5441710885922546
LinearSolver iteration 7, residual = 1.4644985020244865
LinearSolver iteration 8, residual = 1.3521655569577122
LinearSolver iteration 9, residual = 1.336103244954049
LinearSolver iteration 10, residual = 1.259562544285705
LinearSolver iteration 11, residual = 1.1495024170163213
LinearSolver iteration 12, residual = 1.1491221325036693
LinearSolver iteration 13, residual = 1.0583582852278852
LinearSolver iteration 14, residual = 0.9773598464388684
LinearSolver iteration 15, residual = 0.9756045310923616
LinearSolver iteration 16, residual = 0.8670011646769754
LinearSolver iteration 17, residual = 0.7989406664024942
LinearSolver iteration 18, residual = 0.7843579390403251
LinearSolver iteration 19, residual = 0.6959411722935328
LinearSolver iteration 20, residual = 0.6804255085861403
LinearSolver iteration 21, residual = 0.6616381813443352
LinearSolver iteration 22, residual = 0.6248498919661594
LinearSolver iteration 23, residual = 0.6191429972376715
LinearSolver iteration 24, residual = 0.6067983529661145
LinearSolver iteration 25, residual = 0.5794521829418532
LinearSolver iteration 26, residual = 0.5762377937624359
LinearSolver iteration 27, residual = 0.5583795182268424
LinearSolver iteration 28, residual = 0.5344331575565253
LinearSolver iteration 29, residual = 0.5342817943561794
LinearSolver iteration 30, residual = 0.5139643507350377
LinearSolver iteration 31, residual = 0.49698328461047503
LinearSolver iteration 32, residual = 0.4959591001669685
LinearSolver iteration 33, residual = 0.47682299150724106
LinearSolver iteration 34, residual = 0.4730163146750826
LinearSolver iteration 35, residual = 0.4616152701668641
LinearSolver iteration 36, residual = 0.44189221482150465
LinearSolver iteration 37, residual = 0.44184670691877
LinearSolver iteration 38, residual = 0.42122052789831194
LinearSolver iteration 39, residual = 0.39938190425782844
LinearSolver iteration 40, residual = 0.3983788970407686
LinearSolver iteration 41, residual = 0.36309821165343437
LinearSolver iteration 42, residual = 0.3484124801238003
LinearSolver iteration 43, residual = 0.3264000812843651
LinearSolver iteration 44, residual = 0.276430978381519
LinearSolver iteration 45, residual = 0.2764258850664089
LinearSolver iteration 46, residual = 0.22603268049143546
LinearSolver iteration 47, residual = 0.1961534511006646
LinearSolver iteration 48, residual = 0.1863731769161982
LinearSolver iteration 49, residual = 0.13823609788150565
LinearSolver iteration 50, residual = 0.13250126874034732
LinearSolver iteration 51, residual = 0.11243751436651862
LinearSolver iteration 52, residual = 0.0820785306797221
LinearSolver iteration 53, residual = 0.08173368881243895
LinearSolver iteration 54, residual = 0.06039735763812523
LinearSolver iteration 55, residual = 0.044956038457862765
LinearSolver iteration 56, residual = 0.043333358437031726
LinearSolver iteration 57, residual = 0.026912433545661703
LinearSolver iteration 58, residual = 0.023795887364581717
LinearSolver iteration 59, residual = 0.01950604862630428
LinearSolver iteration 60, residual = 0.013335974509596431
LinearSolver iteration 61, residual = 0.013253249776553486
LinearSolver iteration 62, residual = 0.009730686991234623
LinearSolver iteration 63, residual = 0.007615282174440255
LinearSolver iteration 64, residual = 0.007201992058958221
LinearSolver iteration 65, residual = 0.005138225366671139
LinearSolver iteration 66, residual = 0.005061934778723139
LinearSolver iteration 67, residual = 0.004238352789373599
LinearSolver iteration 68, residual = 0.003707802459422031
LinearSolver iteration 69, residual = 0.0036314834557414637
LinearSolver iteration 70, residual = 0.002832986798033658
LinearSolver iteration 71, residual = 0.0026378867951011233
LinearSolver iteration 72, residual = 0.002188349103590497
LinearSolver iteration 73, residual = 0.001522353599047891
LinearSolver iteration 74, residual = 0.0015223533356547355
LinearSolver iteration 75, residual = 0.0010545980209345418
LinearSolver iteration 76, residual = 0.000825252548021231
LinearSolver iteration 77, residual = 0.0007548559927222122
LinearSolver iteration 78, residual = 0.0004631970081715332
LinearSolver iteration 79, residual = 0.0004556031450938567
LinearSolver iteration 80, residual = 0.0003024024909677415
LinearSolver iteration 81, residual = 0.00023793115822170043
LinearSolver iteration 82, residual = 0.0002181038426351963
LinearSolver iteration 83, residual = 0.0001421439069230682
LinearSolver iteration 84, residual = 0.0001311383976718227
LinearSolver iteration 85, residual = 0.00010732819030738014
LinearSolver iteration 86, residual = 6.905427624940651e-05
LinearSolver iteration 87, residual = 6.901033973566872e-05
LinearSolver iteration 88, residual = 4.485409549523471e-05
LinearSolver iteration 89, residual = 3.667425424601136e-05
LinearSolver iteration 90, residual = 3.156119734139944e-05
LinearSolver iteration 91, residual = 2.018654318509908e-05
LinearSolver iteration 92, residual = 1.9643528449685403e-05
LinearSolver iteration 93, residual = 1.3523174721851596e-05
LinearSolver iteration 94, residual = 9.029212123518033e-06
LinearSolver iteration 95, residual = 8.616077389553423e-06
LinearSolver iteration 96, residual = 4.8851277711467706e-06
LinearSolver iteration 97, residual = 4.232314024387145e-06
LinearSolver iteration 98, residual = 3.255348068736631e-06
LinearSolver iteration 99, residual = 1.9291104461704043e-06
LinearSolver iteration 100, residual = 1.924923491320566e-06
LinearSolver iteration 101, residual = 1.1575113412633553e-06
LinearSolver iteration 102, residual = 8.292937123274066e-07
LinearSolver iteration 103, residual = 7.279229052341823e-07
LinearSolver iteration 104, residual = 3.976320693066651e-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.5098463268634976
LinearSolver iteration 3, residual = 2.4713753452091214
LinearSolver iteration 4, residual = 1.6000416541022362
LinearSolver iteration 5, residual = 1.5999941924014407
LinearSolver iteration 6, residual = 1.3318481270718052
LinearSolver iteration 7, residual = 1.2842023741344106
LinearSolver iteration 8, residual = 1.097267216042178
LinearSolver iteration 9, residual = 1.0536946299066943
LinearSolver iteration 10, residual = 0.859140036573202
LinearSolver iteration 11, residual = 0.8019583070152638
LinearSolver iteration 12, residual = 0.6847419608924288
LinearSolver iteration 13, residual = 0.6763718338306801
LinearSolver iteration 14, residual = 0.6146066572025313
LinearSolver iteration 15, residual = 0.6023824596733109
LinearSolver iteration 16, residual = 0.5555606514054204
LinearSolver iteration 17, residual = 0.5430937217136088
LinearSolver iteration 18, residual = 0.5094093630083303
LinearSolver iteration 19, residual = 0.49958044246585603
LinearSolver iteration 20, residual = 0.468723353917065
LinearSolver iteration 21, residual = 0.4561019247305735
LinearSolver iteration 22, residual = 0.42459537484589105
LinearSolver iteration 23, residual = 0.4043814056383923
LinearSolver iteration 24, residual = 0.3455546509233853
LinearSolver iteration 25, residual = 0.3189018314413816
LinearSolver iteration 26, residual = 0.24053411996458152
LinearSolver iteration 27, residual = 0.23388758308230775
LinearSolver iteration 28, residual = 0.13194535898016396
LinearSolver iteration 29, residual = 0.12123099213403957
LinearSolver iteration 30, residual = 0.06382557163320496
LinearSolver iteration 31, residual = 0.06303180309276962
LinearSolver iteration 32, residual = 0.03084994569997178
LinearSolver iteration 33, residual = 0.029346981714241764
LinearSolver iteration 34, residual = 0.015627485603864556
LinearSolver iteration 35, residual = 0.014617906806411023
LinearSolver iteration 36, residual = 0.008292871905729438
LinearSolver iteration 37, residual = 0.008150714306775329
LinearSolver iteration 38, residual = 0.005647994616942128
LinearSolver iteration 39, residual = 0.005620766319648364
LinearSolver iteration 40, residual = 0.003804289007233971
LinearSolver iteration 41, residual = 0.0038014749200121324
LinearSolver iteration 42, residual = 0.002091167783895912
LinearSolver iteration 43, residual = 0.002068155419535695
LinearSolver iteration 44, residual = 0.0010343859929527349
LinearSolver iteration 45, residual = 0.000945083731620937
LinearSolver iteration 46, residual = 0.0005040197865833625
LinearSolver iteration 47, residual = 0.0004947091690997909
LinearSolver iteration 48, residual = 0.0002657759127671414
LinearSolver iteration 49, residual = 0.0002625558560123511
LinearSolver iteration 50, residual = 0.00015169641651776694
LinearSolver iteration 51, residual = 0.00015046497706176342
LinearSolver iteration 52, residual = 9.099101333528002e-05
LinearSolver iteration 53, residual = 9.064467690434781e-05
LinearSolver iteration 54, residual = 5.5494456457668823e-05
LinearSolver iteration 55, residual = 5.442555292178602e-05
LinearSolver iteration 56, residual = 3.374035784882406e-05
LinearSolver iteration 57, residual = 3.31655371468139e-05
LinearSolver iteration 58, residual = 1.8747435764199073e-05
LinearSolver iteration 59, residual = 1.8611134264352855e-05
LinearSolver iteration 60, residual = 9.84861610200471e-06
LinearSolver iteration 61, residual = 9.695503799555885e-06
LinearSolver iteration 62, residual = 5.197401872667092e-06
LinearSolver iteration 63, residual = 5.183534411843179e-06
LinearSolver iteration 64, residual = 2.365090968966256e-06
LinearSolver iteration 65, residual = 2.3016515241079034e-06
LinearSolver iteration 66, residual = 8.885736555415792e-07

Draw (gfu)
Draw (gfp);