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.2774281643545402
LinearSolver iteration 7, residual = 1.2336050111420567
LinearSolver iteration 8, residual = 1.0767167802433257
LinearSolver iteration 9, residual = 0.9764208183923601
LinearSolver iteration 10, residual = 0.8124102954934256
LinearSolver iteration 11, residual = 0.763032628807195
LinearSolver iteration 12, residual = 0.6710519572747343
LinearSolver iteration 13, residual = 0.6468552815219305
LinearSolver iteration 14, residual = 0.596541220502931
LinearSolver iteration 15, residual = 0.5729173657364468
LinearSolver iteration 16, residual = 0.5322237098138078
LinearSolver iteration 17, residual = 0.5244109437623015
LinearSolver iteration 18, residual = 0.49120023781930555
LinearSolver iteration 19, residual = 0.48163025559908507
LinearSolver iteration 20, residual = 0.45495924807777527
LinearSolver iteration 21, residual = 0.4353886755547337
LinearSolver iteration 22, residual = 0.3922054493950524
LinearSolver iteration 23, residual = 0.3702778832189935
LinearSolver iteration 24, residual = 0.301100334830701
LinearSolver iteration 25, residual = 0.2892067537854853
LinearSolver iteration 26, residual = 0.17308383727234025
LinearSolver iteration 27, residual = 0.16928179934913806
LinearSolver iteration 28, residual = 0.08657441731940077
LinearSolver iteration 29, residual = 0.0856289886302004
LinearSolver iteration 30, residual = 0.04685445234257925
LinearSolver iteration 31, residual = 0.046600206745670114
LinearSolver iteration 32, residual = 0.023441789137753227
LinearSolver iteration 33, residual = 0.020118411373025436
LinearSolver iteration 34, residual = 0.013048209565297618
LinearSolver iteration 35, residual = 0.012480171506022018
LinearSolver iteration 36, residual = 0.00913530597535069
LinearSolver iteration 37, residual = 0.008779885649182669
LinearSolver iteration 38, residual = 0.006460732446737279
LinearSolver iteration 39, residual = 0.006058729455791401
LinearSolver iteration 40, residual = 0.004153069889943398
LinearSolver iteration 41, residual = 0.003977244990867421
LinearSolver iteration 42, residual = 0.002549835199294692
LinearSolver iteration 43, residual = 0.002362975425347771
LinearSolver iteration 44, residual = 0.0014897024969474856
LinearSolver iteration 45, residual = 0.0014380081708647685
LinearSolver iteration 46, residual = 0.0009456247234568892
LinearSolver iteration 47, residual = 0.0009042255078697416
LinearSolver iteration 48, residual = 0.0006050248304508075
LinearSolver iteration 49, residual = 0.0005650216216936484
LinearSolver iteration 50, residual = 0.0003967204344210631
LinearSolver iteration 51, residual = 0.0003710361843256435
LinearSolver iteration 52, residual = 0.0002444708028751628
LinearSolver iteration 53, residual = 0.00023962944985237423
LinearSolver iteration 54, residual = 0.00013810579331519193
LinearSolver iteration 55, residual = 0.00012607029920635207
LinearSolver iteration 56, residual = 6.197735441205775e-05
LinearSolver iteration 57, residual = 5.4191891626273075e-05
LinearSolver iteration 58, residual = 2.6817882636778597e-05
LinearSolver iteration 59, residual = 2.3018259029949585e-05
LinearSolver iteration 60, residual = 1.0743945032630732e-05
LinearSolver iteration 61, residual = 1.01986949866314e-05
LinearSolver iteration 62, residual = 4.665057163737831e-06
LinearSolver iteration 63, residual = 4.369820590533748e-06
LinearSolver iteration 64, residual = 1.8770130715748724e-06
LinearSolver iteration 65, residual = 1.524697787081224e-06
LinearSolver iteration 66, residual = 6.16776974871968e-07
LinearSolver iteration 67, residual = 5.436025690540934e-07
LinearSolver iteration 68, residual = 2.509323828701885e-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.0834350054463635
LinearSolver iteration 3, residual = 2.1164861496204943
LinearSolver iteration 4, residual = 1.8366080623996393
LinearSolver iteration 5, residual = 1.5895645511165328
LinearSolver iteration 6, residual = 1.5441710885922568
LinearSolver iteration 7, residual = 1.4644985020244887
LinearSolver iteration 8, residual = 1.352165556957714
LinearSolver iteration 9, residual = 1.3361032449540506
LinearSolver iteration 10, residual = 1.2595625442857061
LinearSolver iteration 11, residual = 1.149502417016322
LinearSolver iteration 12, residual = 1.14912213250367
LinearSolver iteration 13, residual = 1.0583582852278854
LinearSolver iteration 14, residual = 0.9773598464388683
LinearSolver iteration 15, residual = 0.9756045310923614
LinearSolver iteration 16, residual = 0.867001164676975
LinearSolver iteration 17, residual = 0.7989406664024933
LinearSolver iteration 18, residual = 0.7843579390403244
LinearSolver iteration 19, residual = 0.6959411722935317
LinearSolver iteration 20, residual = 0.6804255085861391
LinearSolver iteration 21, residual = 0.661638181344334
LinearSolver iteration 22, residual = 0.6248498919661584
LinearSolver iteration 23, residual = 0.6191429972376705
LinearSolver iteration 24, residual = 0.6067983529661135
LinearSolver iteration 25, residual = 0.5794521829418521
LinearSolver iteration 26, residual = 0.5762377937624348
LinearSolver iteration 27, residual = 0.5583795182268413
LinearSolver iteration 28, residual = 0.5344331575565242
LinearSolver iteration 29, residual = 0.5342817943561783
LinearSolver iteration 30, residual = 0.5139643507350365
LinearSolver iteration 31, residual = 0.49698328461047403
LinearSolver iteration 32, residual = 0.4959591001669675
LinearSolver iteration 33, residual = 0.4768229915072401
LinearSolver iteration 34, residual = 0.47301631467508165
LinearSolver iteration 35, residual = 0.461615270166863
LinearSolver iteration 36, residual = 0.4418922148215033
LinearSolver iteration 37, residual = 0.44184670691876876
LinearSolver iteration 38, residual = 0.4212205278983104
LinearSolver iteration 39, residual = 0.3993819042578266
LinearSolver iteration 40, residual = 0.39837889704076673
LinearSolver iteration 41, residual = 0.3630982116534321
LinearSolver iteration 42, residual = 0.34841248012379805
LinearSolver iteration 43, residual = 0.32640008128436265
LinearSolver iteration 44, residual = 0.2764309783815168
LinearSolver iteration 45, residual = 0.27642588506640675
LinearSolver iteration 46, residual = 0.22603268049144487
LinearSolver iteration 47, residual = 0.1961534511006741
LinearSolver iteration 48, residual = 0.186373176916294
LinearSolver iteration 49, residual = 0.13823609788178992
LinearSolver iteration 50, residual = 0.1325012687403688
LinearSolver iteration 51, residual = 0.11243751436694824
LinearSolver iteration 52, residual = 0.08207853067977816
LinearSolver iteration 53, residual = 0.08173368881247167
LinearSolver iteration 54, residual = 0.06039735763814497
LinearSolver iteration 55, residual = 0.04495603845786216
LinearSolver iteration 56, residual = 0.04333335843703287
LinearSolver iteration 57, residual = 0.026912433545660662
LinearSolver iteration 58, residual = 0.023795887364581245
LinearSolver iteration 59, residual = 0.019506048626301165
LinearSolver iteration 60, residual = 0.01333597450959242
LinearSolver iteration 61, residual = 0.013253249776551372
LinearSolver iteration 62, residual = 0.009730686991208624
LinearSolver iteration 63, residual = 0.007615282174423416
LinearSolver iteration 64, residual = 0.0072019920588607435
LinearSolver iteration 65, residual = 0.005138225366407257
LinearSolver iteration 66, residual = 0.005061934778662059
LinearSolver iteration 67, residual = 0.004238352788097067
LinearSolver iteration 68, residual = 0.003707802458301485
LinearSolver iteration 69, residual = 0.003631483450975063
LinearSolver iteration 70, residual = 0.0028329867651708143
LinearSolver iteration 71, residual = 0.0026378867949023617
LinearSolver iteration 72, residual = 0.0021883488419040713
LinearSolver iteration 73, residual = 0.0015223530715633062
LinearSolver iteration 74, residual = 0.0015223528099513259
LinearSolver iteration 75, residual = 0.0010545922829941416
LinearSolver iteration 76, residual = 0.0008252491396015107
LinearSolver iteration 77, residual = 0.0007548269663237883
LinearSolver iteration 78, residual = 0.00046311468251549647
LinearSolver iteration 79, residual = 0.00045557451462702186
LinearSolver iteration 80, residual = 0.0003019182206067029
LinearSolver iteration 81, residual = 0.00023777684837806537
LinearSolver iteration 82, residual = 0.0002166266953443335
LinearSolver iteration 83, residual = 0.00013837061683604927
LinearSolver iteration 84, residual = 0.00013099270342093993
LinearSolver iteration 85, residual = 9.661669139036816e-05
LinearSolver iteration 86, residual = 6.66615457962145e-05
LinearSolver iteration 87, residual = 6.5880631603458e-05
LinearSolver iteration 88, residual = 4.361716745223726e-05
LinearSolver iteration 89, residual = 3.6661942027680424e-05
LinearSolver iteration 90, residual = 3.1414142104336505e-05
LinearSolver iteration 91, residual = 2.0170567471285503e-05
LinearSolver iteration 92, residual = 1.9638473477387513e-05
LinearSolver iteration 93, residual = 1.3520173836904262e-05
LinearSolver iteration 94, residual = 9.029185359187981e-06
LinearSolver iteration 95, residual = 8.615865572149845e-06
LinearSolver iteration 96, residual = 4.8850951286631905e-06
LinearSolver iteration 97, residual = 4.232312298989344e-06
LinearSolver iteration 98, residual = 3.255344660824827e-06
LinearSolver iteration 99, residual = 1.929110369174175e-06
LinearSolver iteration 100, residual = 1.9249234313914173e-06
LinearSolver iteration 101, residual = 1.1575126033076517e-06
LinearSolver iteration 102, residual = 8.292944938827926e-07
LinearSolver iteration 103, residual = 7.279303693850307e-07
LinearSolver iteration 104, residual = 3.9765585295742756e-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.509846326863497
LinearSolver iteration 3, residual = 2.4713753452091227
LinearSolver iteration 4, residual = 1.6000416541022375
LinearSolver iteration 5, residual = 1.599994192401442
LinearSolver iteration 6, residual = 1.3318481270718097
LinearSolver iteration 7, residual = 1.2842023741344173
LinearSolver iteration 8, residual = 1.097267216042189
LinearSolver iteration 9, residual = 1.0536946299067076
LinearSolver iteration 10, residual = 0.8591400365732205
LinearSolver iteration 11, residual = 0.8019583070152844
LinearSolver iteration 12, residual = 0.6847419608924517
LinearSolver iteration 13, residual = 0.6763718338307039
LinearSolver iteration 14, residual = 0.6146066572025582
LinearSolver iteration 15, residual = 0.60238245967334
LinearSolver iteration 16, residual = 0.5555606514054543
LinearSolver iteration 17, residual = 0.5430937217136457
LinearSolver iteration 18, residual = 0.5094093630083721
LinearSolver iteration 19, residual = 0.4995804424659014
LinearSolver iteration 20, residual = 0.4687233539171173
LinearSolver iteration 21, residual = 0.4561019247306305
LinearSolver iteration 22, residual = 0.42459537484595616
LinearSolver iteration 23, residual = 0.40438140563846403
LinearSolver iteration 24, residual = 0.3455546509234681
LinearSolver iteration 25, residual = 0.31890183144146633
LinearSolver iteration 26, residual = 0.24053411996466448
LinearSolver iteration 27, residual = 0.23388758308238972
LinearSolver iteration 28, residual = 0.13194535898021628
LinearSolver iteration 29, residual = 0.12123099213408998
LinearSolver iteration 30, residual = 0.06382557163323405
LinearSolver iteration 31, residual = 0.06303180309279804
LinearSolver iteration 32, residual = 0.030849945699986447
LinearSolver iteration 33, residual = 0.02934698171425595
LinearSolver iteration 34, residual = 0.015627485603872296
LinearSolver iteration 35, residual = 0.014617906806418115
LinearSolver iteration 36, residual = 0.008292871905733735
LinearSolver iteration 37, residual = 0.008150714306779556
LinearSolver iteration 38, residual = 0.005647994616945142
LinearSolver iteration 39, residual = 0.005620766319651365
LinearSolver iteration 40, residual = 0.0038042890072360245
LinearSolver iteration 41, residual = 0.003801474920014186
LinearSolver iteration 42, residual = 0.0020911677838970706
LinearSolver iteration 43, residual = 0.0020681554195368457
LinearSolver iteration 44, residual = 0.001034385992953319
LinearSolver iteration 45, residual = 0.0009450837316214729
LinearSolver iteration 46, residual = 0.0005040197865836557
LinearSolver iteration 47, residual = 0.0004947091691000772
LinearSolver iteration 48, residual = 0.00026577591276729886
LinearSolver iteration 49, residual = 0.0002625558560125066
LinearSolver iteration 50, residual = 0.00015169641651785907
LinearSolver iteration 51, residual = 0.00015046497706185433
LinearSolver iteration 52, residual = 9.099101333533643e-05
LinearSolver iteration 53, residual = 9.0644676904404e-05
LinearSolver iteration 54, residual = 5.5494456457703674e-05
LinearSolver iteration 55, residual = 5.442555292182028e-05
LinearSolver iteration 56, residual = 3.3740357848845525e-05
LinearSolver iteration 57, residual = 3.316553714683482e-05
LinearSolver iteration 58, residual = 1.8747435764211087e-05
LinearSolver iteration 59, residual = 1.8611134264364687e-05
LinearSolver iteration 60, residual = 9.84861610201142e-06
LinearSolver iteration 61, residual = 9.695503799562402e-06
LinearSolver iteration 62, residual = 5.197401872670926e-06
LinearSolver iteration 63, residual = 5.1835344118470524e-06
LinearSolver iteration 64, residual = 2.3650909689674762e-06
LinearSolver iteration 65, residual = 2.301651524109131e-06
LinearSolver iteration 66, residual = 8.88573655541452e-07

Draw (gfu)
Draw (gfp);