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:

[1]:
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);
[2]:
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:

[3]:
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.

[4]:
(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:

[5]:
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:

[6]:
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:

[7]:
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, mesh)

    return gfu

Higher order Taylor-Hood elements:

[8]:
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:

[9]:
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, mesh)

NgException: UmfpackInverse: Numeric factorization failed.

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

[10]:
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:

[11]:
V = VectorH1(mesh, order=4, dirichlet="wall|inlet|cyl")
# V.SetOrder(TRIG,4)
# V.SetOrder(SEGM,4)
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)
13264
V.ndof = 13264 , Q.ndof = 2406

the mini element:

[12]:
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:

[13]:
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:

[14]:
mp = BilinearForm(p*q*dx).Assemble()

Two right hand sides for the two spaces:

[15]:
f = LinearForm(V).Assemble()
g = LinearForm(Q).Assemble();

Two GridFunctions for velocity and pressure:

[16]:
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.

[17]:
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
[18]:
Draw (gfu);

A stabilised lowest order formulation:

[19]:
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()
[20]:
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
[21]:
Draw (gfu);
[22]:
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.

[23]:
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()
[24]:
# 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.490354257995307
LinearSolver iteration 4, residual = 1.6025072100687436
LinearSolver iteration 5, residual = 1.6024983010872689
LinearSolver iteration 6, residual = 1.33305813757172
LinearSolver iteration 7, residual = 1.28546500251181
LinearSolver iteration 8, residual = 1.0977943618957244
LinearSolver iteration 9, residual = 1.0548387055316295
LinearSolver iteration 10, residual = 0.8601293957401448
LinearSolver iteration 11, residual = 0.8024876442330523
LinearSolver iteration 12, residual = 0.6849389452406784
LinearSolver iteration 13, residual = 0.676610404584322
LinearSolver iteration 14, residual = 0.6147786664099327
LinearSolver iteration 15, residual = 0.602528668180876
LinearSolver iteration 16, residual = 0.5556519286013218
LinearSolver iteration 17, residual = 0.5433107795923555
LinearSolver iteration 18, residual = 0.5095539811129589
LinearSolver iteration 19, residual = 0.4997265271403567
LinearSolver iteration 20, residual = 0.46893652817688825
LinearSolver iteration 21, residual = 0.4562975098173114
LinearSolver iteration 22, residual = 0.4247374804729366
LinearSolver iteration 23, residual = 0.404703864620633
LinearSolver iteration 24, residual = 0.34615234536102096
LinearSolver iteration 25, residual = 0.3190728657926567
LinearSolver iteration 26, residual = 0.24068756246457984
LinearSolver iteration 27, residual = 0.2340879104748613
LinearSolver iteration 28, residual = 0.13210286673685603
LinearSolver iteration 29, residual = 0.12149542221443069
LinearSolver iteration 30, residual = 0.06395535062747418
LinearSolver iteration 31, residual = 0.06313158950830844
LinearSolver iteration 32, residual = 0.030899302683803382
LinearSolver iteration 33, residual = 0.02941571373123929
LinearSolver iteration 34, residual = 0.01569096172127393
LinearSolver iteration 35, residual = 0.01467346862352116
LinearSolver iteration 36, residual = 0.008379750118118865
LinearSolver iteration 37, residual = 0.008238251565404722
LinearSolver iteration 38, residual = 0.005709345997470089
LinearSolver iteration 39, residual = 0.0056849752303122746
LinearSolver iteration 40, residual = 0.0038358687668495085
LinearSolver iteration 41, residual = 0.0038321852291241603
LinearSolver iteration 42, residual = 0.0021014564444746028
LinearSolver iteration 43, residual = 0.0020791718109066885
LinearSolver iteration 44, residual = 0.0010349257957435946
LinearSolver iteration 45, residual = 0.0009452428556792636
LinearSolver iteration 46, residual = 0.0005018308254019071
LinearSolver iteration 47, residual = 0.0004925087878098834
LinearSolver iteration 48, residual = 0.0002624920704936533
LinearSolver iteration 49, residual = 0.00025892211488895417
LinearSolver iteration 50, residual = 0.00014905477022946373
LinearSolver iteration 51, residual = 0.00014768738228125057
LinearSolver iteration 52, residual = 8.918410229413808e-05
LinearSolver iteration 53, residual = 8.884336870670371e-05
LinearSolver iteration 54, residual = 5.443788069477441e-05
LinearSolver iteration 55, residual = 5.330171690493778e-05
LinearSolver iteration 56, residual = 3.3180348230693505e-05
LinearSolver iteration 57, residual = 3.259166338177839e-05
LinearSolver iteration 58, residual = 1.8456984186078528e-05
LinearSolver iteration 59, residual = 1.8306187978107795e-05
LinearSolver iteration 60, residual = 9.729978738601652e-06
LinearSolver iteration 61, residual = 9.554842425629897e-06
LinearSolver iteration 62, residual = 5.1617046055533435e-06
LinearSolver iteration 63, residual = 5.149417860120237e-06
LinearSolver iteration 64, residual = 2.36583603171575e-06
LinearSolver iteration 65, residual = 2.2957609093113194e-06
LinearSolver iteration 66, residual = 8.861914669943986e-07
[25]:
Draw (gfu)
Draw (gfp);
[ ]: