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

69.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.5225660928343492 1.9634170326035583 kappa= 3.7572606786524734
ndof =  217
lammin, lammax= 0.4164879047972151 3.2550708572760314 kappa= 7.815523139527669
ndof =  817
lammin, lammax= 0.37989344758183696 4.308450122441681 kappa= 11.341206724850267
ndof =  3169
lammin, lammax= 0.3560043386452667 5.150172892273618 kappa= 14.466601479835905
ndof =  12481
lammin, lammax= 0.319937819401064 5.825268149129328 kappa= 18.207500945135074
ndof =  49537
lammin, lammax= 0.31667701963457284 6.370873187419749 kappa= 20.117889181764347

ndof =  197377

lammin, lammax= 0.3182030348089053 6.816107452862656 kappa= 21.420623649790215

ndof =  787969

lammin, lammax= 0.3219928450516111 7.183205809357812 kappa= 22.30858827998629

69.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.7651544666507633 0.9987999466057931 kappa= 1.3053572711634862
ndof =  217
lammin, lammax= 0.5892049430232538 0.9967208414315467 kappa= 1.6916369308059416
ndof =  817
lammin, lammax= 0.5497371221971641 0.9970667770354094 kappa= 1.8137155683617994
ndof =  3169
lammin, lammax= 0.4814120552195498 0.996951398840777 kappa= 2.0708899746727645
ndof =  12481
lammin, lammax= 0.45725397386741407 0.9961888045847009 kappa= 2.1786334543121044

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.17521646837299754     
CG iteration 2, residual = 0.00941187124058379     
CG iteration 3, residual = 0.0010818082293433016     
CG iteration 4, residual = 0.0001774124343507029     
CG iteration 5, residual = 3.146247392809796e-05     
CG iteration 6, residual = 5.2014276423085216e-06     
CG iteration 7, residual = 7.973354406340976e-07     
CG iteration 8, residual = 1.6349639615833913e-07     
CG iteration 9, residual = 3.6495382905210624e-08     
CG iteration 10, residual = 6.608662509689858e-09     
CG iteration 11, residual = 1.2577596081448574e-09     
CG iteration 12, residual = 2.468997807675418e-10     
CG iteration 13, residual = 4.94941910837591e-11     
CG iteration 14, residual = 9.056463652129439e-12     
CG iteration 15, residual = 1.7024528929213845e-12     
CG iteration 16, residual = 3.143292625079821e-13     
CG iteration 17, residual = 5.494419632350608e-14     

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

69.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.46135753154810366 0.9962490853348933 kappa= 2.159386196627027

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.17521646837299884     
CG iteration 2, residual = 0.009411871240584515     

CG iteration 3, residual = 0.0010818082293433783     
CG iteration 4, residual = 0.00017741243435071894     
CG iteration 5, residual = 3.146247392809993e-05     
CG iteration 6, residual = 5.20142764230892e-06     
CG iteration 7, residual = 7.973354406342218e-07     
CG iteration 8, residual = 1.6349639615837404e-07     
CG iteration 9, residual = 3.6495382905217427e-08     
CG iteration 10, residual = 6.608662509690732e-09     
CG iteration 11, residual = 1.2577596081450706e-09     
CG iteration 12, residual = 2.46899780767589e-10     
CG iteration 13, residual = 4.9494191083769515e-11     
CG iteration 14, residual = 9.056463652131211e-12     
CG iteration 15, residual = 1.7024528929217936e-12     
CG iteration 16, residual = 3.143292625081005e-13     
CG iteration 17, residual = 5.494419632355397e-14