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.7546 maxew=0.9982 Cond# = 1.323
ndof= 217: minew=0.5814 maxew=0.9966 Cond# = 1.714
ndof= 817: minew=0.5487 maxew=0.9969 Cond# = 1.817
ndof= 3169: minew=0.5072 maxew=0.9970 Cond# = 1.966
ndof= 12481: minew=0.4855 maxew=0.9965 Cond# = 2.052
ndof= 49537: minew=0.4679 maxew=0.9957 Cond# = 2.128
ndof= 197377: minew=0.4648 maxew=0.9949 Cond# = 2.140
ndof= 787969: minew=0.4598 maxew=0.9944 Cond# = 2.163
ndof=3148801: minew=0.4583 maxew=0.9937 Cond# = 2.168
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.17485536352723402
CG iteration 2, residual = 0.010216363008767358
CG iteration 3, residual = 0.0012559902825967464
CG iteration 4, residual = 0.00017268892497020126
CG iteration 5, residual = 2.780650806264886e-05
CG iteration 6, residual = 4.914896534431166e-06
CG iteration 7, residual = 9.226368004868553e-07
CG iteration 8, residual = 1.8365976890034004e-07
CG iteration 9, residual = 3.466556059253677e-08
CG iteration 10, residual = 6.496583860567377e-09
CG iteration 11, residual = 1.1573647765120612e-09
CG iteration 12, residual = 2.602546486822257e-10
CG iteration 13, residual = 5.79547364011061e-11
CG iteration 14, residual = 9.173756057577008e-12
CG iteration 15, residual = 1.7768948746515497e-12
CG iteration 16, residual = 3.7334017096564613e-13
CG iteration 17, residual = 6.852289845642312e-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.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