# 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.

Suppose we have a sequence of hierarchically refined meshes and a sequence of nested finite element spaces

$V_0 \subset V_1 \subset \ldots V_L$

of dimension $$N_l = \operatorname{dim} V_l, l = 0 \ldots L$$. Suppose we also have the prolongation matrices

$P_l \in {\mathbb R}^{N_l \times N_{l-1}}$

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

$A_{l-1} = P_l^T A_l P_l.$

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

$C_{2L}^{-1} = D_l^{-1} + P_l A_{l-1}^{-1} P_l^T.$

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()


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()

[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

$A_{l-1} = P_{l}^T A_l P_l.$
[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()

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