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
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)
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)
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=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:
[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 GridFunction
s 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.4903542579953064
LinearSolver iteration 4, residual = 1.6025072100687432
LinearSolver iteration 5, residual = 1.6024983010872684
LinearSolver iteration 6, residual = 1.3330581375717194
LinearSolver iteration 7, residual = 1.2854650025118097
LinearSolver iteration 8, residual = 1.0977943618957255
LinearSolver iteration 9, residual = 1.054838705531631
LinearSolver iteration 10, residual = 0.8601293957401479
LinearSolver iteration 11, residual = 0.802487644233056
LinearSolver iteration 12, residual = 0.6849389452406832
LinearSolver iteration 13, residual = 0.6766104045843271
LinearSolver iteration 14, residual = 0.614778666409939
LinearSolver iteration 15, residual = 0.6025286681808829
LinearSolver iteration 16, residual = 0.5556519286013301
LinearSolver iteration 17, residual = 0.5433107795923647
LinearSolver iteration 18, residual = 0.5095539811129696
LinearSolver iteration 19, residual = 0.49972652714036836
LinearSolver iteration 20, residual = 0.46893652817690185
LinearSolver iteration 21, residual = 0.45629750981732653
LinearSolver iteration 22, residual = 0.42473748047295384
LinearSolver iteration 23, residual = 0.40470386462065244
LinearSolver iteration 24, residual = 0.3461523453610432
LinearSolver iteration 25, residual = 0.3190728657926797
LinearSolver iteration 26, residual = 0.24068756246460216
LinearSolver iteration 27, residual = 0.23408791047488325
LinearSolver iteration 28, residual = 0.13210286673687255
LinearSolver iteration 29, residual = 0.12149542221444594
LinearSolver iteration 30, residual = 0.06395535062748267
LinearSolver iteration 31, residual = 0.06313158950831695
LinearSolver iteration 32, residual = 0.030899302683807434
LinearSolver iteration 33, residual = 0.02941571373124316
LinearSolver iteration 34, residual = 0.01569096172127595
LinearSolver iteration 35, residual = 0.014673468623523143
LinearSolver iteration 36, residual = 0.00837975011811989
LinearSolver iteration 37, residual = 0.008238251565405736
LinearSolver iteration 38, residual = 0.005709345997470775
LinearSolver iteration 39, residual = 0.00568497523031296
LinearSolver iteration 40, residual = 0.003835868766849986
LinearSolver iteration 41, residual = 0.0038321852291246356
LinearSolver iteration 42, residual = 0.0021014564444748673
LinearSolver iteration 43, residual = 0.0020791718109069513
LinearSolver iteration 44, residual = 0.0010349257957437206
LinearSolver iteration 45, residual = 0.0009452428556793784
LinearSolver iteration 46, residual = 0.0005018308254019656
LinearSolver iteration 47, residual = 0.0004925087878099406
LinearSolver iteration 48, residual = 0.00026249207049368274
LinearSolver iteration 49, residual = 0.00025892211488898334
LinearSolver iteration 50, residual = 0.00014905477022947818
LinearSolver iteration 51, residual = 0.00014768738228126488
LinearSolver iteration 52, residual = 8.918410229414549e-05
LinearSolver iteration 53, residual = 8.884336870671116e-05
LinearSolver iteration 54, residual = 5.4437880694777895e-05
LinearSolver iteration 55, residual = 5.3301716904941116e-05
LinearSolver iteration 56, residual = 3.3180348230694515e-05
LinearSolver iteration 57, residual = 3.259166338177934e-05
LinearSolver iteration 58, residual = 1.8456984186078006e-05
LinearSolver iteration 59, residual = 1.830618797810711e-05
LinearSolver iteration 60, residual = 9.729978738601361e-06
LinearSolver iteration 61, residual = 9.55484242562953e-06
LinearSolver iteration 62, residual = 5.161704605553468e-06
LinearSolver iteration 63, residual = 5.149417860120384e-06
LinearSolver iteration 64, residual = 2.3658360317166303e-06
LinearSolver iteration 65, residual = 2.2957609093121885e-06
LinearSolver iteration 66, residual = 8.861914669951439e-07
[25]:
Draw (gfu)
Draw (gfp);
[ ]: