This page was generated from unit-2.1.5-amg/amg.ipynb.

2.1.5 Algebraic Multigrid Methods

Algebraic multigrid methods (AMG) build a multigrid hierarchy from the given matrix. In contrast to geometric multigrid methods, they do not need a mesh hierarchy. Just one finite element mesh is enough.

AMG takes profit from providing the type of problem (Poisson equation, elasticity, Maxwell, …).

NGSolve comes with builtin AMG solvers for scalar equations, and for Maxwell equations. It provides also interfaces to external, parallel AMG solvers (hypre, gamg, …)

The builtin h1amg

The h1amg preconditioner works for symmetric, scalar problems with nodal degrees of freedom. It uses unsmoothed agglomeration for the generation of coarse spaces.

The coarsening of degrees of freedom is steered by the strength of connections between dofs, one may think of a network of resistors. For this, one finds edge-based weights \(w_E\) such that the energy norm is equivalent to the weighted sum of squared differences:

\[u^T A u \approx \sum_{{\text edges} E} w_E \, (u_{E_1} - u_{E_2})^2\]

\(w_E\) is the edge-weight (the conductivity of each resistor), and \(E_1\) and \(E_2\) are the vertex numbers of the end-points of the edge \(E\). The right hand side is a norm represented by a surrogate matrix \(\tilde A\).

The first task is to determine the edge-weights \(w_E\). If one has access to element-matrices (instead of the assembled matrix), one has better possibilities. One may compute Schur complements with respect to all edges of each element, which gives a surrogate matrix for each element. Then sum up the weights (conductivities) of all elements sharing the edge.

To have access to element matrices, the setup of the surrogate matrix is included into the assembling loop. Thus, the workflow is to

  1. define the biliear-form

  2. define the h1amg preconditioner, which registers at the bilinear-form

  3. finally assemble the bilinear-form, which also runs the setup of the preconditioner

[1]:

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

# minimize memory requirements by switching off tables which we don't need here
import netgen.meshing
netgen.meshing.Mesh.EnableTableClass("edges", False)
netgen.meshing.Mesh.EnableTableClass("faces", False)

with TaskManager():
    mesh = Mesh(unit_cube.GenerateMesh(maxh=0.1))
    for l in range(3): mesh.Refine()

[2]:

# fes = H1(mesh, order=1, order_policy=ORDER_POLICY.CONSTANT)  # todo: fix withtout edge/face tables
fes = FESpace("nodal", mesh, order=1)
print ("ndof=", fes.ndof)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx + 1e-3*u*v*dx)
pre = Preconditioner(a, "h1amg")
with TaskManager():
    a.Assemble();
    lam = EigenValues_Preconditioner(a.mat, pre.mat)
    print (list(lam[0:3]), '...', list(lam[-3:-1]))

ndof= 443525
[0.17997806598214672, 0.20437195578680495, 0.24106066517299723] ... [0.9647917760504281, 0.9927104709152319]

H(curl) - AMG

The hcurlamg is an implementation of the amg from Reitzinger and Schöberl: An algebraic multigrid method for finite element discretizations with edge elements.

It is based on a surrogate matrix for a weighted \(H(\operatorname{curl})\) norm discretized by lowest order Nedelec elements:

\[\| u \|_{L_2, \sigma}^2 + \| \operatorname{curl} u \|_{L_2, \nu}^2 \approx \sum_E w_E \, \Big(\int_E u_{\tau} \Big)^2 + \sum_F w_F \, \Big(\int_F \operatorname{curl}_n u \Big)^2\]

The smoother is a Hiptmair smoother, where a Gauss-Seidel smoother is combined with another Gauss-Seidel smoother for the potential space.

The key is a coarsening which preserves the de Rham sequence over all levels, such that Hiptmair’s smoother is effective also on coarser levels.

More recent, robust coarsening strategies are developed in B. Schwarzenbacher: Robust algebraic solvers for electromagnetics, Master’s Thesis

[3]:

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
# switch on again generation of tables
import netgen.meshing
netgen.meshing.Mesh.EnableTableClass("edges", True)
netgen.meshing.Mesh.EnableTableClass("faces", True)

with TaskManager():
    mesh = Mesh(unit_cube.GenerateMesh(maxh=0.1))
    for l in range(1): mesh.Refine()

[4]:

fes = HCurl(mesh, order=0)
print ("ndof = ", fes.ndof)
u,v = fes.TnT()

a = BilinearForm(curl(u)*curl(v)*dx + 0.01*u*v*dx)

# still crashes from time to time
# pre = Preconditioner(a, "hcurlamg")
# with TaskManager():
#    a.Assemble()
#    lam = EigenValues_Preconditioner(a.mat, pre.mat)
#    print (list(lam[0:3]), '...', list(lam[-3:-1]))

ndof =  57219

[5]:

f = LinearForm(curl(v)[2]*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver

# inv = CGSolver(a.mat, pre.mat, plotrates=False, maxiter=200)
# gfu.vec[:] = inv*f.vec

[ ]: