Preconditioners in NGSolve

Preconditioners are approximative inverses which are used within iterative methods to solve linear or non-linear equations.

Available preconditioners in NGSolve are

  • Jacobi and Block-Jacobi
  • Direct solvers, i.e. sparse factorization
  • Multigrid with different block-smoothers
  • p-version element-level BDDC
In [ ]:
from netgen.csg import unit_cube
from netgen.geom2d import unit_square
from ngsolve import *
import netgen.gui
%gui tk

Define a standard Poisson problem, original mesh-size $h$, polynomial order $p$, refinement levels $l$, and give the name of a preconditioner as argument.

Returns list of tuples of number of unknowns and iterations:

In [ ]:
def SolveProblem(h=0.5, p=1, levels=1, 
                 eliminate_internal = False,
                 precond="local"):
    
    mesh = Mesh(unit_square.GenerateMesh(maxh=h))
    # mesh = Mesh(unit_cube.GenerateMesh(maxh=h))
    fes = H1(mesh, order=p, dirichlet="bottom|left")
    
    u = fes.TrialFunction()
    v = fes.TestFunction()
    a = BilinearForm(fes, flags = { "eliminate_internal" : eliminate_internal })
    a += SymbolicBFI(grad(u)*(grad(v)))
    f = LinearForm(fes)
    f += SymbolicLFI(1*v)
    gfu = GridFunction(fes)
    Draw (gfu)
    c = Preconditioner(a, precond)
    steps = []
    
    for l in range(levels):
        if l > 0: mesh.Refine()
        fes.Update()
        a.Assemble()
        f.Assemble()
        gfu.Update()
        inv = CGSolver(a.mat, c.mat, maxsteps=1000)
        if eliminate_internal:
            f.vec.data += a.harmonic_extension_trans * f.vec
        gfu.vec.data = inv * f.vec
        if eliminate_internal:
            gfu.vec.data += a.harmonic_extension * gfu.vec
            gfu.vec.data += a.inner_solve * f.vec
        steps.append ( (fes.ndof, inv.GetSteps()) )
        Redraw ()
    return steps

The Preconditioner registers itself to the BilinearForm. Whenever the BilinearForm is re-assembled, the Preconditioner is updated as well.

In [ ]:
SolveProblem(precond="local")

The number of CG-iterations with a local (i.e. Jacobi) preconditioner is proportional to $h^{-1} \sim 2^l$:

In [ ]:
res_local = SolveProblem(levels=8, precond="local")
res_local

A geometric multigrid Preconditioner uses the sequence of refined meshes to define a preconditioner of optimal iteration number (and complexity as well). It uses a direct solve on the coarsest level, and block Gauss-Seidel smoothers on the refined levels

In [ ]:
res_mg = SolveProblem(levels=8, precond="multigrid")
res_mg
In [ ]:
import matplotlib.pyplot as plt
plt.xscale("log")
plt.yscale("log")
plt.plot(*zip(*res_local), "-*")
plt.plot(*zip(*res_mg), "-+")
plt.show()

For high order elements we use hierarchical spaces, where the (small) sub-spaces $V_E$, $V_F$, and $V_C$ are generated by basis functions associated with edges, faces, and cells:

$$ V_{hp} = V_{p=1} + \sum_{\text{edges }E} V_E + \sum_{\text{faces }F} V_F + \sum_{\text{cells }C} V_C $$

The system matrix is stored as

$$ A = \left( \begin{array}{cccc} A_{VV} & A_{VE} & A_{VF} & A_{VC} \\ A_{EV} & A_{EE} & A_{EF} & A_{EC} \\ A_{FV} & A_{FE} & A_{FF} & A_{FC} \\ A_{CV} & A_{CE} & A_{CF} & A_{CC} \\ \end{array} \right) $$

The $A_{VV}$-block is exactly the system matrix of a lowest order method.

A multigrid for a high order method uses h-version multigrid for the lowest order block, and only local block-smoothing for the high order bubble basis functions.

In [ ]:
for p in range(1,10):
    r = SolveProblem(h=0.5, p=p, levels=4, eliminate_internal=False, 
                     precond="multigrid")
    print ("p =",p,", res =",r)

We observe that the condition number grows mildly with the order, and remains bounded with mesh refinement. Performing static condensation improves the condition number:

In [ ]:
for p in range(1,10):
    r = SolveProblem(h=0.5, p=p, levels=4, eliminate_internal=True, 
                     precond="multigrid")
    print ("p =",p,", res =",r)

For an element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) one replaces the finite element space by a space connecting only element vertices, but leaving edge and face variables discontinuous. This allows a local elimination of edge and face variables, and thus a cheap global inverse of the replacement matrix. It is used as a preconditioner for the original system matrix:

$$ C_{BDDC}^{-1} = R \widetilde{A}^{-1} R^t $$

Here, $R$ is the averaging operator for the discontinous edge- and face variables.

In contrast to local or multigrid preconditioners, the BDDC - preconditioner needs access to the element matrices. This is exactly the reason why we register the preconditioner with the BilinerForm, and call the bfa.Assemble() later.

In [ ]:
for p in range(1,10):
    r = SolveProblem(h=0.25, p=p, levels=3, eliminate_internal=True, 
                     precond="bddc")
    print ("p =",p,", res =",r)

The BDDC preconditioner needs more iterations, but the work per iteration is less, so performance is similar to multigrid. The BDDC works well in shared memory parallel as well as in distributed memory mode.

In [ ]: