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.4899 maxew=1.9979 Cond# = 4.078
ndof= 217: minew=0.3946 maxew=3.3688 Cond# = 8.537
ndof= 817: minew=0.3824 maxew=4.4660 Cond# = 11.680
ndof= 3169: minew=0.3578 maxew=5.3344 Cond# = 14.909
ndof= 12481: minew=0.3451 maxew=6.0241 Cond# = 17.457
ndof= 49537: minew=0.3479 maxew=6.5755 Cond# = 18.902
ndof= 197377: minew=0.3490 maxew=7.0200 Cond# = 20.113
ndof= 787969: minew=0.3524 maxew=7.3818 Cond# = 20.950
ndof=3148801: minew=0.3549 maxew=7.6793 Cond# = 21.635
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.7671 maxew=0.9981 Cond# = 1.301
ndof= 217: minew=0.5329 maxew=0.9973 Cond# = 1.871
ndof= 817: minew=0.4357 maxew=0.9974 Cond# = 2.289
ndof= 3169: minew=0.3741 maxew=0.9976 Cond# = 2.667
ndof= 12481: minew=0.3201 maxew=0.9972 Cond# = 3.115
ndof= 49537: minew=0.3100 maxew=0.9961 Cond# = 3.213
ndof= 197377: minew=0.2975 maxew=0.9951 Cond# = 3.345
ndof= 787969: minew=0.2922 maxew=0.9943 Cond# = 3.403
ndof=3148801: minew=0.2896 maxew=0.9936 Cond# = 3.431
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.17442004065321884
CG iteration 2, residual = 0.011612915456088663
CG iteration 3, residual = 0.0020605656661244567
CG iteration 4, residual = 0.00033253104993331676
CG iteration 5, residual = 6.718539836727684e-05
CG iteration 6, residual = 1.891733809663533e-05
CG iteration 7, residual = 6.112090250429045e-06
CG iteration 8, residual = 1.914412295208911e-06
CG iteration 9, residual = 5.892236196650654e-07
CG iteration 10, residual = 1.8538787272825326e-07
CG iteration 11, residual = 5.6337454237873476e-08
CG iteration 12, residual = 1.7670272283969393e-08
CG iteration 13, residual = 5.1023284396538875e-09
CG iteration 14, residual = 1.5295254079776071e-09
CG iteration 15, residual = 4.673247097570798e-10
CG iteration 16, residual = 1.43078589469291e-10
CG iteration 17, residual = 4.3865688456266685e-11
CG iteration 18, residual = 1.386520947469764e-11
CG iteration 19, residual = 4.099803119195822e-12
CG iteration 20, residual = 1.1823473229510003e-12
CG iteration 21, residual = 3.5306261714944627e-13
CG iteration 22, residual = 1.1371326103739078e-13
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.4600 maxew=0.9944 Cond# = 2.162
[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.17485526631163548
CG iteration 2, residual = 0.010219366821375229
CG iteration 3, residual = 0.0012591570731032102
CG iteration 4, residual = 0.00017552363474003342
CG iteration 5, residual = 2.97080869626909e-05
CG iteration 6, residual = 5.680101091137623e-06
CG iteration 7, residual = 1.0923124855353862e-06
CG iteration 8, residual = 2.1670664862334042e-07
CG iteration 9, residual = 3.80424301803236e-08
CG iteration 10, residual = 6.4597819081562295e-09
CG iteration 11, residual = 1.3067246266885869e-09
CG iteration 12, residual = 3.1816356237853314e-10
CG iteration 13, residual = 5.2925528874192093e-11
CG iteration 14, residual = 9.784534280154563e-12
CG iteration 15, residual = 2.1453264228410207e-12
CG iteration 16, residual = 3.8312527256316257e-13
CG iteration 17, residual = 7.308497426731703e-14