This page was generated from unit-2.1.3-multigrid/multigrid.ipynb.
2.1.3 Multigrid and Multilevel Methods¶
In the previous tutorial \(\S\)2.1.1, we saw how to construct and use a geometric multigrid preconditioner using a keyword argument to the Preconditioner. This tutorial delves deeper into the construction of such preconditioners. You will see how additive and multiplicative preconditioners can be implemented efficiently by recursive algorithms using NGSolve’s python interface.
Additive Multilevel Preconditioner¶
Suppose we have a sequence of hierarchically refined meshes and a sequence of nested finite element spaces
of dimension \(N_l = \operatorname{dim} V_l, l = 0 \ldots L\). Suppose we also have the prolongation matrices
with the property that \(\underline v_l = P_l \underline v_{l-1}\) where \(\underline v_k\) is the vector of coefficients in the finite element basis expansion of a \(v_k \in V_k\). For the stiffness matrix \(A_k\) of a Galerkin discretization, there holds
Let \(D_l = \operatorname{diag} A_l\) be the Jacobi preconditioner (or some similar, cheap and local preconditioner).
The construction of the preconditioner is motivated by the theory of 2-level preconditioners. The 2-level preconditioner involving levels \(l-1\) and level \(l\), namely
has optimal condition number by the additive Schwarz lemma. However, a direct inversion of the matrix \(A_{l-1}\) (is up to constant factor) as expensive as the inversion of \(A_l\).
The ideal of the multilevel preconditioner is to replace that inverse by a recursion: \begin{eqnarray*} C_{ML,0}^{-1} & = & A_0^{-1} \\ C_{ML,l}^{-1} & = & D_l^{-1} + P_l C_{ML,l-1}^{-1} P_l^T \qquad \text{for} \; l = 1, \ldots L \end{eqnarray*}
We can implement this using a python class with a recursive constructor. We show two implementation techniques, one in long form, and one in short form.
[1]:
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
This is the long form implementation, which is a useful starting point for improvising to more complex preconditioners you might need for other purposes.
[2]:
class MLPreconditioner(BaseMatrix):
def __init__ (self, fes, level, mat, coarsepre):
super().__init__()
self.fes = fes
self.level = level
self.mat = mat
self.coarsepre = coarsepre
if level > 0:
self.localpre = mat.CreateSmoother(fes.FreeDofs())
else:
self.localpre = mat.Inverse(fes.FreeDofs())
def Mult (self, x, y):
"""Return y = C[l] * x = (inv(D) + P * C[l-1] * P.T) * x
when l>0 and inv(A[0]) * x when l=0. """
if self.level == 0:
y.data = self.localpre * x
return
hx = x.CreateVector(copy=True)
self.fes.Prolongation().Restrict(self.level, hx) # hx <- P.T * x
cdofs = self.fes.Prolongation().LevelDofs(self.level-1)
y[cdofs] = self.coarsepre * hx[cdofs] # y = C[l-1] * hx
self.fes.Prolongation().Prolongate(self.level, y) # y <- P * C[l-1] * P.T * x
y += self.localpre * x # y += inv(D) * x
def Shape (self):
return self.localpre.shape
def CreateVector (self, col):
return self.localpre.CreateVector(col)
Here is the considerably shorter version implementing the same preconditioner using operator notation.
[3]:
def MLPreconditioner2(fes, level, mat, coarsepre):
prol = fes.Prolongation().Operator(level) # get P
localpre = mat.CreateSmoother(fes.FreeDofs()) # get inv(D)
return localpre + prol @ coarsepre @ prol.T # Return: inv(D) + P * C[l-1] * P.T * x
[4]:
a.Assemble() # Make exact inverse at current/coarsest level
pre = a.mat.Inverse(fes.FreeDofs())
[5]:
for l in range(9):
mesh.Refine()
a.Assemble()
pre = MLPreconditioner(fes,l+1, a.mat, pre)
lam = EigenValues_Preconditioner(a.mat, pre)
print("ndof=%7d: minew=%.4f maxew=%1.4f Cond# = %5.3f"
%(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))
ndof= 61: minew=0.4873 maxew=1.9985 Cond# = 4.102
ndof= 217: minew=0.3931 maxew=3.3696 Cond# = 8.572
ndof= 817: minew=0.3817 maxew=4.4677 Cond# = 11.703
ndof= 3169: minew=0.3525 maxew=5.3370 Cond# = 15.140
ndof= 12481: minew=0.3381 maxew=6.0275 Cond# = 17.826
ndof= 49537: minew=0.3418 maxew=6.5797 Cond# = 19.251
ndof= 197377: minew=0.3483 maxew=7.0252 Cond# = 20.172
ndof= 787969: minew=0.3531 maxew=7.3879 Cond# = 20.920
ndof=3148801: minew=0.3569 maxew=7.6863 Cond# = 21.535
Multiplicative Multigrid Preconditioning¶
The multigrid preconditioner combines the local preconditioner and the coarse grid step multiplicatively or sequentially. The multigrid preconditioning action is defined as follows:
Multigrid \(C_{MG,l}^{-1} : d \mapsto w\)
if l = 0, set \(w = A_0^{-1} d\) and return
w = 0
presmoothing, \(m_l\) steps:
\[w \leftarrow w + D_{pre}^{-1} (d - A w)\]coasre grid correction:
\[w \leftarrow w + P_l C_{MG,l-1}^{-1} P_l^{T} (d - A w)\]postsmoothing, \(m_l\) steps:
\[w \leftarrow w + D_{post}^{-1} (d - A w)\]
If the preconditioners \(D_{pre}\) and \(D_{post}\) from the pre-smoothing and post-smoothing iterations are transposed to each other, the overall preconditioner is symmetric. If the pre- and post-smoothing iterations are non-expansive, the overall preconditioner is positive definite. These properties can be realized using forward Gauss-Seidel for pre-smoothing, and backward Gauss-Seidel for post-smoothing. Here is its implementation as a recursive python class.
[6]:
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
[7]:
class MGPreconditioner(BaseMatrix):
def __init__ (self, fes, level, mat, coarsepre):
super().__init__()
self.fes = fes
self.level = level
self.mat = mat
self.coarsepre = coarsepre
if level > 0:
self.localpre = mat.CreateSmoother(fes.FreeDofs())
else:
self.localpre = mat.Inverse(fes.FreeDofs())
def Mult (self, d, w):
if self.level == 0:
w.data = self.localpre * d
return
prol = self.fes.Prolongation().Operator(self.level)
w[:] = 0
self.localpre.Smooth(w,d)
res = d - self.mat * w
w += prol @ self.coarsepre @ prol.T * res
self.localpre.SmoothBack(w,d)
def Shape (self):
return self.localpre.shape
def CreateVector (self, col):
return self.localpre.CreateVector(col)
[8]:
a.Assemble()
pre = MGPreconditioner(fes, 0, a.mat, None)
for l in range(9):
mesh.Refine()
a.Assemble()
pre = MGPreconditioner(fes,l+1, a.mat, pre)
lam = EigenValues_Preconditioner(a.mat, pre)
print("ndof=%7d: minew=%.4f maxew=%1.4f Cond# = %5.3f"
%(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))
ndof= 61: minew=0.7538 maxew=0.9982 Cond# = 1.324
ndof= 217: minew=0.5839 maxew=0.9967 Cond# = 1.707
ndof= 817: minew=0.5580 maxew=0.9969 Cond# = 1.787
ndof= 3169: minew=0.5155 maxew=0.9969 Cond# = 1.934
ndof= 12481: minew=0.4819 maxew=0.9964 Cond# = 2.068
ndof= 49537: minew=0.4624 maxew=0.9957 Cond# = 2.153
ndof= 197377: minew=0.4600 maxew=0.9951 Cond# = 2.163
ndof= 787969: minew=0.4552 maxew=0.9945 Cond# = 2.185
ndof=3148801: minew=0.4535 maxew=0.9939 Cond# = 2.192
A quick conjugate gradient call on this finest grid gives you an idea of typical number of iterations and practial efficiency.
[9]:
f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec
CG iteration 1, residual = 0.17479885997159333
CG iteration 2, residual = 0.010318399831683133
CG iteration 3, residual = 0.00119629114962409
CG iteration 4, residual = 0.00015784111089870974
CG iteration 5, residual = 2.436840504910884e-05
CG iteration 6, residual = 4.218743956723545e-06
CG iteration 7, residual = 7.844578876896069e-07
CG iteration 8, residual = 1.5730309835747064e-07
CG iteration 9, residual = 3.05125711973268e-08
CG iteration 10, residual = 6.050151991063636e-09
CG iteration 11, residual = 1.2603156795373275e-09
CG iteration 12, residual = 2.882358568550667e-10
CG iteration 13, residual = 4.942863170316655e-11
CG iteration 14, residual = 7.930127327478935e-12
CG iteration 15, residual = 1.8311213548968712e-12
CG iteration 16, residual = 3.5775265113237374e-13
CG iteration 17, residual = 6.570052414101423e-14
Projection matrices from the finest level¶
It is often not feasible to assemble matrices on the coarse level, for example when solving non-linear problems. Then it is useful to calculate coarse grid matrices from the matrix on the finest level using the Galerkin property
[10]:
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
for l in range(8):
mesh.Refine()
a.Assemble(); # only matrix at the finest level provided to mg
[11]:
class ProjectedMG(BaseMatrix):
def __init__ (self, fes, mat, level):
super(ProjectedMG, self).__init__()
self.fes = fes
self.level = level
self.mat = mat
if level > 0:
self.prol = fes.Prolongation().CreateMatrix(level)
self.rest = self.prol.CreateTranspose()
coarsemat = self.rest @ mat @ self.prol # multiply matrices
self.localpre = mat.CreateSmoother(fes.FreeDofs())
self.coarsepre = ProjectedMG(fes, coarsemat, level-1)
else:
self.localpre = mat.Inverse(fes.FreeDofs())
def Mult (self, d, w):
if self.level == 0:
w.data = self.localpre * d
return
w[:] = 0
self.localpre.Smooth(w,d)
res = d - self.mat * w
w += self.prol @ self.coarsepre @ self.rest * res
self.localpre.SmoothBack(w,d)
def Shape (self):
return self.localpre.shape
def CreateVector (self, col):
return self.localpre.CreateVector(col)
[12]:
pre = ProjectedMG(fes, a.mat, fes.mesh.levels-1)
[13]:
lam = EigenValues_Preconditioner(a.mat, pre)
print("ndof=%7d: minew=%.4f maxew=%1.4f Cond# = %5.3f"
%(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))
ndof= 787969: minew=0.4555 maxew=0.9945 Cond# = 2.183
[14]:
f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec
CG iteration 1, residual = 0.17479876359909335
CG iteration 2, residual = 0.010321542152665342
CG iteration 3, residual = 0.0011995832870883441
CG iteration 4, residual = 0.000160704993156155
CG iteration 5, residual = 2.62115695170076e-05
CG iteration 6, residual = 4.904951915285568e-06
CG iteration 7, residual = 9.312891708592661e-07
CG iteration 8, residual = 1.8398650620547436e-07
CG iteration 9, residual = 3.336911330328518e-08
CG iteration 10, residual = 6.665725793884852e-09
CG iteration 11, residual = 1.576559823212486e-09
CG iteration 12, residual = 3.0441306818647564e-10
CG iteration 13, residual = 4.730217015638341e-11
CG iteration 14, residual = 1.0015891413412979e-11
CG iteration 15, residual = 2.0799708154904935e-12
CG iteration 16, residual = 3.8274874773413833e-13
CG iteration 17, residual = 7.975171227542562e-14