1. Is NGSolve getting Lazy?#
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
1.1. Some people like the SolvePDE function#
A left-over from pre-Python NGSolve times, but useful …
Code from i-tutorials 1.3, going back to UM3
fes = H1(mesh, order=2, dirichlet="left|right") # Dirichlet boundaries in space
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
gfu.Set(sin(y), BND) # interpolate Dirichlet boundary conditions
c = preconditioners.Local(a)
c.Update() # may or may not need update
Takes care of Dirichlet boundary conditions, static condensation, …
solvers.BVP(bf=a, lf=f, gf=gfu, pre=c)
Draw(gfu);
CG iteration 1, residual = 1.7445591828021156
CG iteration 2, residual = 0.9453692996140045
CG iteration 3, residual = 0.821243134657245
CG iteration 4, residual = 0.5740477176991614
CG iteration 5, residual = 0.5334023460549775
CG iteration 6, residual = 0.4010588114326537
CG iteration 7, residual = 0.45157986816168266
CG iteration 8, residual = 0.4976855831892793
CG iteration 9, residual = 0.41942622092149906
CG iteration 10, residual = 0.3884084339952495
CG iteration 11, residual = 0.2565169299453453
CG iteration 12, residual = 0.15370841654363404
CG iteration 13, residual = 0.09767549832929286
CG iteration 14, residual = 0.06778307402838821
CG iteration 15, residual = 0.054832410959067816
CG iteration 16, residual = 0.047967454244897845
CG iteration 17, residual = 0.052956025122294914
CG iteration 18, residual = 0.04851085763073881
CG iteration 19, residual = 0.027217241881460782
CG iteration 20, residual = 0.01845451815778772
CG iteration 21, residual = 0.013002049663921851
CG iteration 22, residual = 0.008618425151477321
CG iteration 23, residual = 0.006724423125332083
CG iteration 24, residual = 0.005350075748960324
CG iteration 25, residual = 0.0038539739612640177
CG iteration 26, residual = 0.002472361881551127
CG iteration 27, residual = 0.0017281235463580912
CG iteration 28, residual = 0.0012737305828211995
CG iteration 29, residual = 0.0008434892203301497
CG iteration 30, residual = 0.00052712070205152
CG iteration 31, residual = 0.00030847495405016834
CG iteration 32, residual = 0.00022121819588806217
CG iteration 33, residual = 0.00016649771319683117
CG iteration 34, residual = 0.00011279274220861203
CG iteration 35, residual = 8.022756901757145e-05
CG iteration 36, residual = 5.53128700988198e-05
CG iteration 37, residual = 3.334353398465425e-05
CG iteration 38, residual = 2.03858473072738e-05
CG iteration 39, residual = 1.2840021780916396e-05
CG iteration 40, residual = 8.836700317527927e-06
CG iteration 41, residual = 5.510608306599427e-06
CG iteration 42, residual = 3.4402195972713406e-06
CG iteration 43, residual = 2.288579471937665e-06
CG iteration 44, residual = 1.5202450527317397e-06
CG iteration 45, residual = 1.1133006213701104e-06
CG iteration 46, residual = 7.533479990323252e-07
CG iteration 47, residual = 4.921358808507492e-07
CG iteration 48, residual = 3.448624502415105e-07
CG iteration 49, residual = 2.16571736999024e-07
CG iteration 50, residual = 1.446500617374522e-07
CG iteration 51, residual = 9.460134405177877e-08
CG iteration 52, residual = 6.527780382544902e-08
CG iteration 53, residual = 4.0282206754630344e-08
CG iteration 54, residual = 2.164922275560863e-08
CG iteration 55, residual = 1.2752669980955774e-08
1.2. Version 2#
from UM5:
More verbose defining linear or non-linear problems:
fes = H1(mesh, order=2, dirichlet="left|right")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
c = preconditioners.Local(a)
f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
Solve(a*gfu==f, pre=c, \
dirichlet=mesh.BoundaryCF({'left' : sin(y), 'right':-sin(y)}))
Draw (gfu);
1.3. New attempt#
Only creator objects go into the Solve function
u,v = H1(mesh, order=2).TnT()
lhs = grad(u)*grad(v)*dx
rhs = 1*v*dx
gfu = Solve(lhs==rhs, u[BND("left")]==sin(y), u[BND("right")]==-sin(y), \
preconditioners.Local(), solvers.CGSolver.Creator())
Draw (gfu);
If we want to give arguments, provide it to the creator:
gfu = Solve(lhs==rhs, u[BND("left|right")]==sin(y), \
preconditioners.Local(blocktype="vertexpatch"),
solvers.CGSolver.Creator(maxiter=40, printrates=True),
verbose=2)
Draw (gfu);
PreconditionerCreator: <class 'ngsolve.comp.LocalPreconditioner'>
Dirichlet: 0: 0101
CG iteration 1, residual = 1.7733437901620588
CG iteration 2, residual = 1.1453273468733678
CG iteration 3, residual = 0.6866019245445
CG iteration 4, residual = 0.7060033921272163
CG iteration 5, residual = 0.47609401985320576
CG iteration 6, residual = 0.4609895460437994
CG iteration 7, residual = 0.42955072878998857
CG iteration 8, residual = 0.4838173447672771
CG iteration 9, residual = 0.6225022179287477
CG iteration 10, residual = 0.35212321873310704
CG iteration 11, residual = 0.1763226013931562
CG iteration 12, residual = 0.09575249786255077
CG iteration 13, residual = 0.10073374322633914
CG iteration 14, residual = 0.06365792411386856
CG iteration 15, residual = 0.06701818823812081
CG iteration 16, residual = 0.06405096932036038
CG iteration 17, residual = 0.05167905933457604
CG iteration 18, residual = 0.0412980009661957
CG iteration 19, residual = 0.024886524598506245
CG iteration 20, residual = 0.019083497886155532
CG iteration 21, residual = 0.012807776159985722
CG iteration 22, residual = 0.011159423036993006
CG iteration 23, residual = 0.007757893483453212
CG iteration 24, residual = 0.005915255718843835
CG iteration 25, residual = 0.0036259891514376793
CG iteration 26, residual = 0.002886547195973667
CG iteration 27, residual = 0.0016569332826158595
CG iteration 28, residual = 0.0010948614742171967
CG iteration 29, residual = 0.000758775759483668
CG iteration 30, residual = 0.0005398160033887696
CG iteration 31, residual = 0.0003904589435242948
CG iteration 32, residual = 0.00030263064484697163
CG iteration 33, residual = 0.00021730187023039684
CG iteration 34, residual = 0.00010743880346134202
CG iteration 35, residual = 7.537532022605711e-05
CG iteration 36, residual = 3.547994186365287e-05
CG iteration 37, residual = 2.582523266269845e-05
CG iteration 38, residual = 1.735201236460509e-05
CG iteration 39, residual = 1.3386074712006034e-05
CG iteration 40, residual = 7.4808381275877675e-06
WARNING: CG did not converge to TOL
1.4. What are the new components?#
1.4.1. Region descriptors#
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=1)
u,v = fes.TnT()
gfu = GridFunction(fes)
gfu[BND("left")] = y
gfu[BND("top")] = 1-x
Draw (gfu);
type (BND("left"))
ngsolve.comp.RegionDescriptor
print (BND("left"))
BND(left)
print (u[BND("left")])
<ngsolve.comp.DirichletBoundary object at 0x10d7701b0>
print (u[BND("left")]==3)
<ngsolve.comp.DirichletBC object at 0x10d770bf0>
reg = mesh[BND("left")]
print ("reg =", reg)
print ("reg.Mask() =", reg.Mask())
reg = <ngsolve.comp.Region object at 0x10d771070>
reg.Mask() = 0: 0001
1.4.2. Our three versions of setting boundary values of a GridFunction#
The Set function performs element-wise \(L_2\)-projection, and then averaging coupling dofs.
It only needs basis-functions.
For historic reasons, the Set first zeros the complete GridFunction, thus consecutive calls to Set overwrites previous values.
gfu = GridFunction(fes)
vals = mesh.BoundaryCF( { 'left' : y**2, 'right' : -y } )
gfu.Set(vals, definedon=mesh.Boundaries("left|right"))
Draw (gfu);
As we got degrees-of-freedom into NGSolve (aka DualShapes), we could implement canonical interpolation. This is now available for most spaces. The Interpolate function sets values on the given region, in leaves degrees of freedom not in the region unchanged. Note that the maximal value is now precisely as expected.
gfu = GridFunction(fes)
gfu.Interpolate(y**2, definedon=mesh.Boundaries("right"))
gfu.Interpolate(-y, definedon=mesh.Boundaries("left"))
Draw (gfu);
The mesh enters via the GridFunction, and via the Region as well. There is some redundancy.
The new syntax using GridFunction.__setitem__is a compact notation for Interpolate:
gfu = GridFunction(fes)
gfu[BND("left")] = -y
gfu[BND("right")] = y
Draw (gfu);
1.5. Objects and Object creators#
If we have a bilinear-from, we create a preconditioner, with some flags
a = BilinearForm(grad(u)*grad(v)*dx)
c = preconditioners.Local(a, GS=True)
type(c)
ngsolve.comp.LocalPreconditioner
print(c)
type = Local Preconditioner
But if we don’t have a BilinearForm object yet, we cannot create a preconditioner.
First attempt: Give a Creator static function to the preconditioner class, it provides a factory/builder/creator class
creator = preconditioners.MultiGrid.Creator(cycle=2, smoother="block")
type(creator)
ngsolve.comp.PreconditionerCreator
pre = creator(a)
type(pre), print (pre)
Multigrid preconditioner
bilinear-form = biform_from_py
smoothertype = block
(ngsolve.comp.MultiGridPreconditioner, None)
Second attempt: Use the same class as class creator and final class object:
pre_creator = preconditioners.MultiGrid(cycle=2, smoother="block")
type (pre_creator)
ngsolve.comp.MultiGridPreconditioner
We did not give a BilinearForm object, so we have an incomplete object. However, it knows what type of (creator-)object it is, and stores arguments
pre_creator.IsCreator()
True
finally, we can use it to create the real preconditioner object:
pre = pre_creator.Create(a, cycle=2)
type (pre)
ngsolve.comp.MultiGridPreconditioner
pre.IsCreator()
False
The code of the new Solve-function (WIP):
class VariationalEquationSolver:
def __init__(self, equation, *args : DirichletBC | Preconditioner, verbose=0, **kwargs):
self.bf = BilinearForm(equation.igls)
self.dirichlet = [a for a in args if isinstance(a, DirichletBC)]
self.fes = self.bf.space
self.mesh = self.fes.mesh
if self.dirichlet:
self.dreg = sum((self.mesh[d.vbn] for d in self.dirichlet),
self.mesh[self.dirichlet[0].vbn])
for a in args:
if (verbose>=5): print ("got argument of type", type(a))
if isinstance(a, Preconditioner):
if a.IsCreator():
if verbose>=2: print ("PreconditionerCreator: ", type(a))
if self.dirichlet:
if verbose>=2: print ("Dirichlet: ", self.dreg.Mask())
self.pre = a.Create(self.bf, additional_dirichlet_constraints=self.dreg)
else:
self.pre = a.Create(self.bf)
else:
self.pre = a
if isinstance(a, LinearSolverCreator):
self.linear_solver_creator = a
if isinstance(a, SparseFactorizationCreator):
if verbose>=2: print ("SparseFactorizationCreator: ", type(a))
self.sparse_factorization_creator = a
if pre := kwargs.get('pre'):
self.pre=pre
def Solve(self):
with TaskManager():
gf = GridFunction(self.fes)
for d in self.dirichlet:
gf.ComponentFromProxy(d.proxy)[d.vbn] = d.val
self.bf.AssembleLinearization(gf.vec)
if hasattr(self, 'linear_solver_creator'):
inv = self.linear_solver_creator(self.bf.mat, self.pre)
else:
freedofs = self.fes.FreeDofs()
if self.dirichlet:
freedofs = freedofs&(~self.fes.GetDofs(self.dreg))
if hasattr(self, 'sparse_factorization_creator'):
inv = self.sparse_factorization_creator(self.bf.mat, freedofs)
else:
inv = self.bf.mat.Inverse(freedofs)
gf.vec.data -= inv * self.bf.Apply(gf.vec)
return gf
def SolveVE (equation, *args, **kwargs):
return VariationalEquationSolver(equation,*args,**kwargs).Solve()