88. H(curl) - AMG

We present 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 edge-weights stem from the \(L_{2,\sigma}\) norm. One could take the diagonal of the edge-element mass matrix, or also Schur complements with respect to the edges. The \(\operatorname{curl}\)-semi-norm is represented by the face terms. The weights can be computed by lumping the lowest-order Raviart-Thomas mass matrix.

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.

\[\begin{split} \!\!\!\!\!\!\!\!\! \begin{array}{ccccccc} & B_{\operatorname{grad}} & & B_{\operatorname{curl}} & & B_{\operatorname{div}} & \\[-0.5em] V^v & -\!\!\!\longrightarrow & V^e & -\!\!\!-\!\!\!\longrightarrow & V^f & -\!\!\!-\!\!\!\longrightarrow & V^c \\[1em] \;\;\;\;\downarrow \Pi^W & & \;\;\;\;\downarrow \Pi^V & & \;\;\;\;\downarrow \Pi^Q & & \;\;\;\;\downarrow \Pi^S \\[1em] & B_{\operatorname{grad}} & & B_{\operatorname{curl}} & & B_{\operatorname{div}} & \\[-0.5em] V^v_{coarse} & -\!\!\!\longrightarrow & V^e_{coarse} & -\!\!\!-\!\!\!\longrightarrow & V^f_{coarse} & -\!\!\!-\!\!\!\longrightarrow & V^c_{coarse} \end{array} \end{split}\]

The coarsening of edges is derived from coarsening of vertices. \(E_{IJ}\) is a coarse grid edge if and only if \(I \neq J\), and there are fine grid vertices \(i\) and \(j\) s.t.:

\[ I = Ind(i), \quad J = Ind(j), \qquad E_{ij} \mbox{ is a fine grid edge} \]

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

89. H(curl) - AMG in NGSolve

This amg method is implemented as hcurlamg preconditioner in NGSolve.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw

coil = Cylinder( Axes( (0,0,-0.4), Z, X), h=0.8,r=0.4) \
    - Cylinder( Axes( (0,0,-0.4), Z, X), h=0.8,r=0.2)
box = Box( (-2,-2,-2), (2,2,2) )

coil.faces.col=(1,0,0)
coil.faces.maxh=0.1
coil.solids.name="coil"
box.faces.col=(0,0,1,0.3)
box.faces.name="outer"
air = box - coil
shape = Glue( [coil,air] )
Draw (shape);

mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.2)) # .Curve(3)
Draw (mesh);

fes = HCurl(mesh, order=2, nograds=True) # , dirichlet="outer")
print ("ndof=", fes.ndof)
u,v = fes.TnT()

with TaskManager():
    mu = 4*pi*1e-7
    a = BilinearForm(1/mu*curl(u)*curl(v)*dx + 1e-6/mu*u*v*dx \
                     + 1e8*u.Trace()*v.Trace()*ds("outer"))
    # pre = preconditioners.HCurlAMG(a)
    pre = preconditioners.MultiGrid(a, coarsetype="hcurlamg", coarseflags={ "steps" : 5 }) 
    a.Assemble()
    f = LinearForm( CF((y,-x,0))*v*dx("coil", bonus_intorder=5)).Assemble()

ndof= 212102

WARNING: kwarg 'coarseflags' is an undocumented flags option for class <class 'ngsolve.comp.MultiGridPreconditioner'>, maybe there is a typo?

gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver

with TaskManager():
    inv = CGSolver(a.mat, pre.mat, plotrates=True, maxiter=100)
    gfu.vec[:] = inv*f.vec

<Figure size 640x480 with 0 Axes>

ea = { "euler_angles" : (-60, 0, -11) }
clipping = { "clipping" : { "y":1, "z":0, "dist":0.5} }

s = 0.1*2
N = 10
p = [(-s+2*s*i/N,-s+2*s*j/N,-s+2*s*k/N) for i in range(1,N) for j in range(1,N) for k in range(1,N)]

fieldlines = curl(gfu)._BuildFieldLines(mesh, p, num_fieldlines=N**3//5, randomized=True, length=2)

Draw(curl(gfu), mesh,  "X", draw_vol=False, draw_surf=True, objects=[fieldlines], \
     min=0, max=1e-8, autoscale=False, settings={"Objects": {"Surface": False}},
    **ea, **clipping);

help (preconditioners.MultiGrid)

Help on class MultiGridPreconditioner in module ngsolve.comp:

class MultiGridPreconditioner(Preconditioner)
 |   Keyword arguments can be:
 |  inverse:
 |    Inverse type used in Preconditioner.
 |  test: bool = False
 |    Computes condition number for preconditioner, if testout file
 |    is set, prints eigenvalues to file.
 |  block:
 |    use block-Jacobi/block-Gauss-Seidel
 |  blocktype: str = vertexpatch
 |    Blocktype used in compound FESpace for smoothing
 |    blocks. Options: vertexpatch, edgepatch
 |  updateall: bool = False
 |    Update all smoothing levels when calling Update
 |  smoother: string = 'point'
 |    Smoother between multigrid levels, available options are:
 |      'point': Gauss-Seidel-Smoother
 |      'line':  Anisotropic smoother
 |      'block': Block smoother
 |  coarsetype: string = direct
 |    How to solve coarse problem.
 |  coarsesmoothingsteps: int = 1
 |    If coarsetype is smoothing, then how many smoothingsteps will be done.
 |  updatealways: bool = False
 |
 |  Method resolution order:
 |      MultiGridPreconditioner
 |      Preconditioner
 |      ngsolve.la.BaseMatrix
 |      NGS_Object
 |      pybind11_builtins.pybind11_object
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  SetDirectSolverCluster(...)
 |      SetDirectSolverCluster(self: ngsolve.comp.MultiGridPreconditioner, arg0: list) -> None
 |
 |  __init__(...)
 |      __init__(self: ngsolve.comp.MultiGridPreconditioner, bf: ngsolve.comp.BilinearForm, name: str = 'multigrid', lo_preconditioner: Optional[ngsolve.comp.Preconditioner] = None, **kwargs) -> None
 |
 |  ----------------------------------------------------------------------
 |  Static methods defined here:
 |
 |  __flags_doc__(...) method of builtins.PyCapsule instance
 |      __flags_doc__() -> dict
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |
 |  __dict__
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from Preconditioner:
 |
 |  Test(...)
 |      Test(self: ngsolve.comp.Preconditioner) -> None
 |
 |  Update(...)
 |      Update(self: ngsolve.comp.Preconditioner) -> None
 |
 |      Update preconditioner
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from Preconditioner:
 |
 |  mat
 |      matrix of the preconditioner
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from ngsolve.la.BaseMatrix:
 |
 |  AsVector(...)
 |      AsVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
 |
 |      Interprets the matrix values as a vector
 |
 |  CreateColVector(...)
 |      CreateColVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
 |
 |  CreateDeviceMatrix(...)
 |      CreateDeviceMatrix(self: ngsolve.la.BaseMatrix) -> BaseMatrix
 |
 |  CreateMatrix(...)
 |      CreateMatrix(self: ngsolve.la.BaseMatrix) -> BaseMatrix
 |
 |      Create matrix of same dimension and same sparsestructure
 |
 |  CreateRowVector(...)
 |      CreateRowVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
 |
 |  CreateVector(...)
 |      CreateVector(self: ngsolve.la.BaseMatrix, colvector: bool = False) -> ngsolve.la.BaseVector
 |
 |  GetInverseType(...)
 |      GetInverseType(self: ngsolve.la.BaseMatrix) -> str
 |
 |  GetOperatorInfo(...)
 |      GetOperatorInfo(self: ngsolve.la.BaseMatrix) -> str
 |
 |  Inverse(...)
 |      Inverse(self: ngsolve.la.BaseMatrix, freedofs: pyngcore.pyngcore.BitArray = None, inverse: Optional[str] = None, flags: pyngcore.pyngcore.Flags = <pyngcore.pyngcore.Flags object at 0x1089645f0>) -> BaseMatrix
 |
 |      Calculate inverse of sparse matrix
 |      Parameters:
 |
 |      freedofs : BitArray
 |        If set, invert only the rows/columns the matrix defined by the bit array, otherwise invert the whole matrix
 |
 |      inverse : string
 |        Solver to use, allowed values are:
 |          sparsecholesky - internal solver of NGSolve for symmetric matrices
 |          umfpack        - solver by Suitesparse/UMFPACK (if NGSolve was configured with USE_UMFPACK=ON)
 |          pardiso        - PARDISO, either provided by libpardiso (USE_PARDISO=ON) or Intel MKL (USE_MKL=ON).
 |                           If neither Pardiso nor Intel MKL was linked at compile-time, NGSolve will look
 |                           for libmkl_rt in LD_LIBRARY_PATH (Unix) or PATH (Windows) at run-time.
 |
 |  Mult(...)
 |      Mult(self: ngsolve.la.BaseMatrix, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |  MultAdd(...)
 |      MultAdd(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. MultAdd(self: ngsolve.la.BaseMatrix, value: float, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |      2. MultAdd(self: ngsolve.la.BaseMatrix, value: complex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |  MultScale(...)
 |      MultScale(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. MultScale(self: ngsolve.la.BaseMatrix, value: float, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |      2. MultScale(self: ngsolve.la.BaseMatrix, value: complex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |  MultTrans(...)
 |      MultTrans(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. MultTrans(self: ngsolve.la.BaseMatrix, value: float, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |      2. MultTrans(self: ngsolve.la.BaseMatrix, value: complex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |  MultTransAdd(...)
 |      MultTransAdd(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. MultTransAdd(self: ngsolve.la.BaseMatrix, value: float, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |      2. MultTransAdd(self: ngsolve.la.BaseMatrix, value: complex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
 |
 |  ToDense(...)
 |      ToDense(self: ngsolve.la.BaseMatrix) -> object
 |
 |  __add__(...)
 |      __add__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
 |
 |  __iadd__(...)
 |      __iadd__(self: ngsolve.la.BaseMatrix, mat: ngsolve.la.BaseMatrix) -> None
 |
 |  __matmul__(...)
 |      __matmul__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
 |
 |  __mul__(...)
 |      __mul__(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. __mul__(self: BaseMatrix, arg0: ngsolve.la.BaseVector) -> ngla::DynamicVectorExpression
 |
 |      2. __mul__(self: BaseMatrix, arg0: ngsolve.la.MultiVector) -> ngsolve.la.MultiVectorExpr
 |
 |  __neg__(...)
 |      __neg__(self: BaseMatrix) -> BaseMatrix
 |
 |  __radd__(...)
 |      __radd__(self: BaseMatrix, arg0: int) -> BaseMatrix
 |
 |  __rmul__(...)
 |      __rmul__(*args, **kwargs)
 |      Overloaded function.
 |
 |      1. __rmul__(self: BaseMatrix, value: float) -> BaseMatrix
 |
 |      2. __rmul__(self: BaseMatrix, value: complex) -> BaseMatrix
 |
 |  __str__(...)
 |      __str__(self: ngsolve.la.BaseMatrix) -> str
 |
 |  __sub__(...)
 |      __sub__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
 |
 |  __timing__(...)
 |      __timing__(self: ngsolve.la.BaseMatrix, runs: int = 10) -> float
 |
 |  matvec(...)
 |      matvec(self: BaseMatrix, arg0: ngsolve.la.BaseVector) -> ngsolve.la.BaseVector
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from ngsolve.la.BaseMatrix:
 |
 |  H
 |      Return conjugate transpose of matrix (WIP, only partially supported)
 |
 |  T
 |      Return transpose of matrix
 |
 |  comm
 |
 |  dtype
 |
 |  height
 |      Height of the matrix
 |
 |  is_complex
 |      is the matrix complex-valued ?
 |
 |  local_mat
 |
 |  nze
 |      number of non-zero elements
 |
 |  shape
 |
 |  width
 |      Width of the matrix
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from NGS_Object:
 |
 |  __memory__
 |
 |  flags
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from NGS_Object:
 |
 |  name
 |
 |  ----------------------------------------------------------------------
 |  Static methods inherited from pybind11_builtins.pybind11_object:
 |
 |  __new__(*args, **kwargs) class method of pybind11_builtins.pybind11_object
 |      Create and return a new object.  See help(type) for accurate signature.