85. Multigrid and Multilevel Methods

Multigrid (MG) and Multilevel (ML) algorithms provide preconditioners with optimal condition numbers \(\kappa (C^{-1} A) = O(1)\), and optimal computational complexity \(O(N)\).

They can be seen as extension of the two level overlapping domain decomposition method to more levels.

• Ulrich Trottenberg, Cornelius W. Oosterlee, Anton Schuller: Multigrid, Academic Press, 2001

• Wolfgang Hackbusch: Multi-Grid Methods and Applications, Springer, 1985

In short, both methods are sub-space correction methods where the space splitting is defined by all basis functions together from a hierarchy of grids. Multilevel methods are additive Schwarz methods, while Multigrid methods are the multiplicative Schwarz methods. Both can be implemented efficiently by recursive algorithms.

We will present different theories for their analysis, one is based on sub-space decomposition using the ASM lemma, the other one uses the smoothing-and-approximation properties.

85.1. Multilevel preconditioner

We have a sequence of hierarchically refined meshes, which lead to 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\). We think of mesh-sizes \(h_l = 2^{-l}\), leading to finite element space dimensions \(2^{dl}\), where \(d\) is the spatial dimension.

We have prolongation matrices

\[ P_l \in {\mathbb R}^{N_l \times N_{l-1}}. \]

If \(v_{l-1}\) is a finite element function in \(V_{l-1}\) represented by the coefficient vector \(\underline v_{l-1}\). Then \(\underline v_l = P_l \underline v_{l-1}\) is the coefficient vector of the same function represented by the basis functions of \(V_l\).

If \(A_l\) and \(A_{l-1}\) are discretization matrices by a Galerkin method, then 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).

2-level method

A 2-level preconditioner involving levels \(l-1\) and level \(l\) is

\[ C_{2L}^{-1} = D_l^{-1} + P_l A_{l-1}^{-1} P_l^T \]

By the ASM - Lemma we have the tools to analyze that this preconditioner has optimal condition number. However, a direct inversion of the matrix \(A_{l-1}\) is, up to a constant factor, as expensive as the inversion of \(A_l\). The multilevel method is to replace the inverses recursively by a multilevel preconditioner. On the coarsest level we invert the matrix \(A_0\):

\[\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*}\]

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)

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):
        if self.level == 0:
            y.data = self.localpre * x
            return
        hx = x.CreateVector(copy=True)
        self.fes.Prolongation().Restrict(self.level, hx)
        cdofs = self.fes.Prolongation().LevelDofs(self.level-1)
        y[cdofs] = self.coarsepre * hx[cdofs] 
        self.fes.Prolongation().Prolongate(self.level, y)
        y += self.localpre * x

    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

With operator notation we can define the multilevel preconditioner also as follows:

def MLPreconditioner2(fes, level, mat, coarsepre):
    prol = fes.Prolongation().Operator(level)
    localpre = mat.CreateSmoother(fes.FreeDofs())
    return localpre + prol @ coarsepre @ prol.T

a.Assemble()
pre = a.mat.Inverse(fes.FreeDofs())

for l in range(8):
    mesh.Refine()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MLPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.5041881887985011 1.9554176933071035 kappa= 3.8783488720093486
ndof =  217
lammin, lammax= 0.4077070793994908 3.2339129373496 kappa= 7.931951885929501
ndof =  817
lammin, lammax= 0.37039417252180995 4.278495350080584 kappa= 11.551195098320973
ndof =  3169
lammin, lammax= 0.3445249894360645 5.1149130463560315 kappa= 14.846275896353333
ndof =  12481
lammin, lammax= 0.3253496279478505 5.787442674704374 kappa= 17.78837956941518
ndof =  49537
lammin, lammax= 0.3221458682366649 6.333364925485618 kappa= 19.65992908787768

ndof =  197377
lammin, lammax= 0.3225863521811482 6.782386381082096 kappa= 21.025025811611062

ndof =  787969

lammin, lammax= 0.3260334941433799 7.156685776004728 kappa= 21.950768569984497

85.2. Multigrid Preconditioning

The multigrid preconditioner combines the local preconditioner and the coarse grid step sequentially. The multigrid preconditioning action is define 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:

  \[\hat w = w + D_{pre}^{-1} (d - A w)\]

• coasre grid correction:

  \[ \hat w = w + P_l C_{MG,l-1}^{-1} P_l^{T} (d - A w) \]

• postsmoothing, \(m_l\) steps:

  \[\hat w = 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. This can be obtained by applying forward Gauss-Seidel for pre-smoothing, and backward Gauss-Seidel for post-smoothing.

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

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):
        # print ("level = ", self.level)
        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)

a.Assemble()
pre = MGPreconditioner(fes, 0, a.mat, None)

for l in range(5):
    mesh.Refine()
    fes.Update()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MGPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.7670141189280819 0.9985799776429505 kappa= 1.3019056012143384
ndof =  217
lammin, lammax= 0.5982873742488635 0.9966894497582455 kappa= 1.665904200317727
ndof =  817
lammin, lammax= 0.5751097730546433 0.9971514409408411 kappa= 1.733845411189175
ndof =  3169
lammin, lammax= 0.4942050326345251 0.9969910302502516 kappa= 2.017363167945615
ndof =  12481
lammin, lammax= 0.4589480199019186 0.9962363630079554 kappa= 2.1706954160535656

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.175135158473104     
CG iteration 2, residual = 0.009653901012112616     
CG iteration 3, residual = 0.001056372809227722     
CG iteration 4, residual = 0.0001718886616474045     
CG iteration 5, residual = 3.0817777740443086e-05     
CG iteration 6, residual = 5.198477096536569e-06     
CG iteration 7, residual = 9.044940575146323e-07     
CG iteration 8, residual = 2.0336348749728586e-07     
CG iteration 9, residual = 4.124150542645893e-08     
CG iteration 10, residual = 7.144541562621555e-09     
CG iteration 11, residual = 1.3230901428300202e-09     
CG iteration 12, residual = 2.6990549474525664e-10     
CG iteration 13, residual = 5.5716845089930815e-11     
CG iteration 14, residual = 1.0692928806682936e-11     
CG iteration 15, residual = 2.0932080370146976e-12     
CG iteration 16, residual = 4.2863020673126905e-13     
CG iteration 17, residual = 8.849795358713877e-14     

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

85.3. 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 \]

[ Requires NGSolve updates from 14.04.2021, still tricky: freedofs on coarser levels]

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(5):
     mesh.Refine()
    
a.Assemble();

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)

# print (fes.mesh.levels)
pre = ProjectedMG(fes, a.mat, fes.mesh.levels-1)

print ("ndof: ", fes.ndof)
lam = EigenValues_Preconditioner(a.mat, pre)
print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof:  12481
lammin, lammax= 0.46089246840624576 0.9962665529200481 kappa= 2.1616030228593495

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

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

CG iteration 1, residual = 0.17513515847310573     
CG iteration 2, residual = 0.009653901012113338     
CG iteration 3, residual = 0.0010563728092277838     
CG iteration 4, residual = 0.00017188866164740763     
CG iteration 5, residual = 3.081777774044404e-05     
CG iteration 6, residual = 5.1984770965365126e-06     
CG iteration 7, residual = 9.044940575145592e-07     
CG iteration 8, residual = 2.0336348749726746e-07     
CG iteration 9, residual = 4.1241505426457495e-08     
CG iteration 10, residual = 7.144541562621548e-09     
CG iteration 11, residual = 1.3230901428299404e-09     
CG iteration 12, residual = 2.699054947452313e-10     
CG iteration 13, residual = 5.5716845089927913e-11     
CG iteration 14, residual = 1.069292880668287e-11     
CG iteration 15, residual = 2.093208037014634e-12     
CG iteration 16, residual = 4.2863020673123654e-13     
CG iteration 17, residual = 8.849795358713399e-14