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 = 3424 , Q.ndof = 2406

---------------------------------------------------------------------------
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 = 9840 , Q.ndof = 7218

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

5028
V.ndof = 5028 , Q.ndof = 2406

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.518937179332287
LinearSolver iteration 2, residual = 2.3020782066183405
LinearSolver iteration 3, residual = 1.9171806179628295
LinearSolver iteration 4, residual = 1.5318309921095468
LinearSolver iteration 5, residual = 1.4827182650518633
LinearSolver iteration 6, residual = 1.2760882423711855
LinearSolver iteration 7, residual = 1.2316328218317878
LinearSolver iteration 8, residual = 1.0758256973668074
LinearSolver iteration 9, residual = 0.9742346042735229
LinearSolver iteration 10, residual = 0.8116311768045531
LinearSolver iteration 11, residual = 0.7620621123363509
LinearSolver iteration 12, residual = 0.6711057838761307
LinearSolver iteration 13, residual = 0.6467734480369526
LinearSolver iteration 14, residual = 0.5963885550482779
LinearSolver iteration 15, residual = 0.5715760562984827
LinearSolver iteration 16, residual = 0.5314172872757961
LinearSolver iteration 17, residual = 0.523188572014509
LinearSolver iteration 18, residual = 0.48916660118256133
LinearSolver iteration 19, residual = 0.48026612300248067
LinearSolver iteration 20, residual = 0.45410105844800014
LinearSolver iteration 21, residual = 0.4339175906792935
LinearSolver iteration 22, residual = 0.38970115024248825
LinearSolver iteration 23, residual = 0.3665030310531088
LinearSolver iteration 24, residual = 0.29515807082938256
LinearSolver iteration 25, residual = 0.2840989870359845
LinearSolver iteration 26, residual = 0.17263506434502293
LinearSolver iteration 27, residual = 0.16533278871221124
LinearSolver iteration 28, residual = 0.09207002076023608
LinearSolver iteration 29, residual = 0.08806383521916354
LinearSolver iteration 30, residual = 0.045707773020785346
LinearSolver iteration 31, residual = 0.04558872526001847
LinearSolver iteration 32, residual = 0.022509806937679445
LinearSolver iteration 33, residual = 0.01911902981935344
LinearSolver iteration 34, residual = 0.012572671188669397
LinearSolver iteration 35, residual = 0.012232920688064624
LinearSolver iteration 36, residual = 0.008796244991863617
LinearSolver iteration 37, residual = 0.008540624451742073
LinearSolver iteration 38, residual = 0.006168751068326351
LinearSolver iteration 39, residual = 0.005859614741848614
LinearSolver iteration 40, residual = 0.003919721191004725
LinearSolver iteration 41, residual = 0.0037309231359978446
LinearSolver iteration 42, residual = 0.002409644869021594
LinearSolver iteration 43, residual = 0.002233246013776026
LinearSolver iteration 44, residual = 0.001440117277538793
LinearSolver iteration 45, residual = 0.0014169553919431872
LinearSolver iteration 46, residual = 0.0009214338668519829
LinearSolver iteration 47, residual = 0.0008998464054884886
LinearSolver iteration 48, residual = 0.0005762503566586648
LinearSolver iteration 49, residual = 0.0005472865054142066
LinearSolver iteration 50, residual = 0.0003758901019645419
LinearSolver iteration 51, residual = 0.0003524418202292021
LinearSolver iteration 52, residual = 0.00023429185631029594
LinearSolver iteration 53, residual = 0.000231577683020344
LinearSolver iteration 54, residual = 0.00013261675061463022
LinearSolver iteration 55, residual = 0.00012528252842354204
LinearSolver iteration 56, residual = 5.852579935273627e-05
LinearSolver iteration 57, residual = 5.220921278375006e-05
LinearSolver iteration 58, residual = 2.445981833622924e-05
LinearSolver iteration 59, residual = 2.1011457947562368e-05
LinearSolver iteration 60, residual = 9.894725207855162e-06
LinearSolver iteration 61, residual = 9.325347229145796e-06
LinearSolver iteration 62, residual = 4.39262521268901e-06
LinearSolver iteration 63, residual = 4.195388653320418e-06
LinearSolver iteration 64, residual = 1.7610904076427593e-06
LinearSolver iteration 65, residual = 1.487707702456541e-06
LinearSolver iteration 66, residual = 5.821123043092815e-07
LinearSolver iteration 67, residual = 5.153650968031574e-07
LinearSolver iteration 68, residual = 2.3534539744561125e-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.103717015785737
LinearSolver iteration 2, residual = 3.1240529736545715
LinearSolver iteration 3, residual = 2.132679123744692
LinearSolver iteration 4, residual = 1.8523744559707935
LinearSolver iteration 5, residual = 1.5949503379856076
LinearSolver iteration 6, residual = 1.5480493182937696
LinearSolver iteration 7, residual = 1.4639703507160722
LinearSolver iteration 8, residual = 1.3511179026398876
LinearSolver iteration 9, residual = 1.3380334785486305
LinearSolver iteration 10, residual = 1.2598116080932027
LinearSolver iteration 11, residual = 1.1487182058160144
LinearSolver iteration 12, residual = 1.1487081802610593
LinearSolver iteration 13, residual = 1.0553572664456368
LinearSolver iteration 14, residual = 0.976593791636269
LinearSolver iteration 15, residual = 0.9718806996194856
LinearSolver iteration 16, residual = 0.8599661205136232
LinearSolver iteration 17, residual = 0.8005914699770541
LinearSolver iteration 18, residual = 0.7815317337264195
LinearSolver iteration 19, residual = 0.695091815685787
LinearSolver iteration 20, residual = 0.681097261054522
LinearSolver iteration 21, residual = 0.6601652993538863
LinearSolver iteration 22, residual = 0.6238986254835396
LinearSolver iteration 23, residual = 0.6195450703831556
LinearSolver iteration 24, residual = 0.6055276228835791
LinearSolver iteration 25, residual = 0.5784318059307156
LinearSolver iteration 26, residual = 0.5763221221836843
LinearSolver iteration 27, residual = 0.5571394201839083
LinearSolver iteration 28, residual = 0.5339897335990003
LinearSolver iteration 29, residual = 0.5339897224407189
LinearSolver iteration 30, residual = 0.5122111603632727
LinearSolver iteration 31, residual = 0.49709386368401864
LinearSolver iteration 32, residual = 0.49522266954245237
LinearSolver iteration 33, residual = 0.4757510640129083
LinearSolver iteration 34, residual = 0.4731826388330148
LinearSolver iteration 35, residual = 0.4601388105401166
LinearSolver iteration 36, residual = 0.4412886936788078
LinearSolver iteration 37, residual = 0.44120267904223176
LinearSolver iteration 38, residual = 0.4187428770452127
LinearSolver iteration 39, residual = 0.39925869154183435
LinearSolver iteration 40, residual = 0.3965165497660975
LinearSolver iteration 41, residual = 0.35958302389869107
LinearSolver iteration 42, residual = 0.34918705384313015
LinearSolver iteration 43, residual = 0.32081334730538613
LinearSolver iteration 44, residual = 0.27496623675114984
LinearSolver iteration 45, residual = 0.2739395154602775
LinearSolver iteration 46, residual = 0.22065410237349553
LinearSolver iteration 47, residual = 0.19723595462027937
LinearSolver iteration 48, residual = 0.18105301698811066
LinearSolver iteration 49, residual = 0.13692381837465906
LinearSolver iteration 50, residual = 0.13287558442797537
LinearSolver iteration 51, residual = 0.10926674275131437
LinearSolver iteration 52, residual = 0.08098202942409063
LinearSolver iteration 53, residual = 0.08094489870875758
LinearSolver iteration 54, residual = 0.056389394203058316
LinearSolver iteration 55, residual = 0.04514989779565494
LinearSolver iteration 56, residual = 0.040324817474304074
LinearSolver iteration 57, residual = 0.02515404088463173
LinearSolver iteration 58, residual = 0.023960521236286676
LinearSolver iteration 59, residual = 0.018012559701725786
LinearSolver iteration 60, residual = 0.01300583626679384
LinearSolver iteration 61, residual = 0.012973772374450137
LinearSolver iteration 62, residual = 0.009169395609592765
LinearSolver iteration 63, residual = 0.00763667146865482
LinearSolver iteration 64, residual = 0.0068485873124133475
LinearSolver iteration 65, residual = 0.004943988286513263
LinearSolver iteration 66, residual = 0.004911909351469885
LinearSolver iteration 67, residual = 0.004081919199534889
LinearSolver iteration 68, residual = 0.0035790619585567927
LinearSolver iteration 69, residual = 0.0034957389244193965
LinearSolver iteration 70, residual = 0.002702516402196035
LinearSolver iteration 71, residual = 0.0026133921922323797
LinearSolver iteration 72, residual = 0.002044199670776769
LinearSolver iteration 73, residual = 0.0014952781332316647
LinearSolver iteration 74, residual = 0.0014711480663878792
LinearSolver iteration 75, residual = 0.000991617872172705
LinearSolver iteration 76, residual = 0.0008714858820776782
LinearSolver iteration 77, residual = 0.0007056267834200474
LinearSolver iteration 78, residual = 0.0004773645624879543
LinearSolver iteration 79, residual = 0.00046424789573374227
LinearSolver iteration 80, residual = 0.00028094430610380993
LinearSolver iteration 81, residual = 0.0002548453786863786
LinearSolver iteration 82, residual = 0.00020256676252789458
LinearSolver iteration 83, residual = 0.00013267945882768027
LinearSolver iteration 84, residual = 0.00013194521868198242
LinearSolver iteration 85, residual = 8.593924383152689e-05
LinearSolver iteration 86, residual = 6.239009989375969e-05
LinearSolver iteration 87, residual = 5.9815516703039264e-05
LinearSolver iteration 88, residual = 3.849464222210089e-05
LinearSolver iteration 89, residual = 3.2605363484935703e-05
LinearSolver iteration 90, residual = 2.8322245608267862e-05
LinearSolver iteration 91, residual = 1.82690807881192e-05
LinearSolver iteration 92, residual = 1.8051591698283513e-05
LinearSolver iteration 93, residual = 1.1906749700372199e-05
LinearSolver iteration 94, residual = 8.918371402129329e-06
LinearSolver iteration 95, residual = 7.708397821999984e-06
LinearSolver iteration 96, residual = 4.360196894806923e-06
LinearSolver iteration 97, residual = 4.241381799516438e-06
LinearSolver iteration 98, residual = 2.799363495931647e-06
LinearSolver iteration 99, residual = 1.8482147124704675e-06
LinearSolver iteration 100, residual = 1.7777878936095626e-06
LinearSolver iteration 101, residual = 9.742813141308085e-07
LinearSolver iteration 102, residual = 8.494431817617426e-07
LinearSolver iteration 103, residual = 6.2071818814675e-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.21762787865306
LinearSolver iteration 2, residual = 3.5470429131109387
LinearSolver iteration 3, residual = 2.4903542579953064
LinearSolver iteration 4, residual = 1.6025072100687436
LinearSolver iteration 5, residual = 1.6024983010872689
LinearSolver iteration 6, residual = 1.3330581375717196
LinearSolver iteration 7, residual = 1.2854650025118095
LinearSolver iteration 8, residual = 1.0977943618957242
LinearSolver iteration 9, residual = 1.0548387055316286
LinearSolver iteration 10, residual = 0.8601293957401428
LinearSolver iteration 11, residual = 0.8024876442330491
LinearSolver iteration 12, residual = 0.6849389452406744
LinearSolver iteration 13, residual = 0.6766104045843175
LinearSolver iteration 14, residual = 0.614778666409927
LinearSolver iteration 15, residual = 0.6025286681808695
LinearSolver iteration 16, residual = 0.5556519286013136
LinearSolver iteration 17, residual = 0.5433107795923466
LinearSolver iteration 18, residual = 0.5095539811129488
LinearSolver iteration 19, residual = 0.4997265271403462
LinearSolver iteration 20, residual = 0.46893652817687664
LinearSolver iteration 21, residual = 0.45629750981729916
LinearSolver iteration 22, residual = 0.4247374804729232
LinearSolver iteration 23, residual = 0.404703864620619
LinearSolver iteration 24, residual = 0.3461523453610062
LinearSolver iteration 25, residual = 0.31907286579264144
LinearSolver iteration 26, residual = 0.2406875624645667
LinearSolver iteration 27, residual = 0.23408791047484778
LinearSolver iteration 28, residual = 0.1321028667368506
LinearSolver iteration 29, residual = 0.12149542221442378
LinearSolver iteration 30, residual = 0.06395535062747042
LinearSolver iteration 31, residual = 0.06313158950830462
LinearSolver iteration 32, residual = 0.03089930268380123
LinearSolver iteration 33, residual = 0.029415713731237348
LinearSolver iteration 34, residual = 0.01569096172127279
LinearSolver iteration 35, residual = 0.014673468623520265
LinearSolver iteration 36, residual = 0.008379750118118187
LinearSolver iteration 37, residual = 0.00823825156540408
LinearSolver iteration 38, residual = 0.005709345997469587
LinearSolver iteration 39, residual = 0.005684975230311783
LinearSolver iteration 40, residual = 0.0038358687668491494
LinearSolver iteration 41, residual = 0.0038321852291237977
LinearSolver iteration 42, residual = 0.0021014564444743703
LinearSolver iteration 43, residual = 0.002079171810906462
LinearSolver iteration 44, residual = 0.0010349257957434621
LinearSolver iteration 45, residual = 0.0009452428556791472
LinearSolver iteration 46, residual = 0.0005018308254018377
LinearSolver iteration 47, residual = 0.0004925087878098165
LinearSolver iteration 48, residual = 0.000262492070493613
LinearSolver iteration 49, residual = 0.00025892211488891465
LinearSolver iteration 50, residual = 0.00014905477022943823
LinearSolver iteration 51, residual = 0.0001476873822812253
LinearSolver iteration 52, residual = 8.918410229412129e-05
LinearSolver iteration 53, residual = 8.884336870668705e-05
LinearSolver iteration 54, residual = 5.443788069476345e-05
LinearSolver iteration 55, residual = 5.330171690492707e-05
LinearSolver iteration 56, residual = 3.318034823068621e-05
LinearSolver iteration 57, residual = 3.259166338177118e-05
LinearSolver iteration 58, residual = 1.845698418607426e-05
LinearSolver iteration 59, residual = 1.8306187978103492e-05
LinearSolver iteration 60, residual = 9.729978738599421e-06
LinearSolver iteration 61, residual = 9.554842425627775e-06
LinearSolver iteration 62, residual = 5.161704605552153e-06
LinearSolver iteration 63, residual = 5.149417860119061e-06
LinearSolver iteration 64, residual = 2.36583603171532e-06
LinearSolver iteration 65, residual = 2.295760909310915e-06
LinearSolver iteration 66, residual = 8.861914669940874e-07

Draw (gfu)
Draw (gfp);