Vertex patch smoothing for p-multigrid preconditioners for the Poisson problem

In this tutorial, we explore using PETSc PC as a building block inside the NGSolve preconditioning infrastructure. Not all the preconditioning strategies are equally effective for the discretisations of the Poisson problem here considered. This demo is intended to provide a starting point for the exploration of preconditioning strategies rather than providing a definitive answer. We begin by creating a discretisation of the Poisson problem using H1 elements, in particular, we consider the usual variational formulation

\[\text{find } u\in H^1_0(\Omega) \text{ s.t. } a(u,v) := \int_{\Omega} \nabla u\cdot \nabla v \; d\vec{x} = L(v) := \int_{\Omega} fv\; d\vec{x}\qquad v\in H^1_0(\Omega).\]

Such a discretisation can easily be constructed using NGSolve as follows:

from ngsolve import *
import netgen.gui
import netgen.meshing as ngm
from mpi4py.MPI import COMM_WORLD

if COMM_WORLD.rank == 0:
   ngmesh = unit_square.GenerateMesh(maxh=0.1)
   for _ in range(4):
      ngmesh.Refine()
   mesh = Mesh(ngmesh.Distribute(COMM_WORLD))
else:
   mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
order = 2
fes = H1(mesh, order=order, dirichlet="left|right|top|bottom")
print("Number of degrees of freedom: ", fes.ndof)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
f = LinearForm(fes)
f += 32 * (y*(1-y)+x*(1-x)) * v * dx
a.Assemble()
f.Assemble()

We now construct an NGSolve preconditioner wrapping a PETSc PC object. In particular, we begin considering an algebraic multigrid preconditioner constructed using HYPRE. We will use the preconditioner just discussed inside NGSolve’s own linear algebra solvers.

from ngsPETSc.pc import *
from ngsolve.krylovspace import CG
pre = Preconditioner(a, "PETScPC", pc_type="hypre")
gfu = GridFunction(fes)
print("-------------------|HYPRE p={}|-------------------".format(order))
gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pre.mat, printrates=True)
Draw(gfu)

We see that the HYPRE preconditioner is quite effective for the Poisson problem discretised using linear elements, but it is not as effective for higher-order elements.

Preconditioners performance HYPRE

Preconditioner

p=1

p=2

p=3

p=4

HYPRE

10 (1.57e-12)

70 (2.22e-13)

42 (3.96e-13)

70 (1.96e-13)

To overcome this issue we will use a two-level additive Schwarz preconditioner. In this case, we will use as fine space correction the inverse of the stiffness matrices associated with the patch of a vertex.

def VertexPatchBlocks(mesh, fes):
 blocks = []
 freedofs = fes.FreeDofs()
 for v in mesh.vertices:
     vdofs = set()
     for el in mesh[v].elements:
         vdofs |= set(d for d in fes.GetDofNrs(el)
                      if freedofs[d])
     blocks.append(vdofs)
 return blocks

blocks = VertexPatchBlocks(mesh, fes)
dofs = BitArray(fes.ndof); dofs[:] = True
blockjac = ASMPreconditioner(a.mat, dofs, blocks=blocks,
                               solverParameters={"pc_type": "asm",
                                                 "sub_ksp_type": "preonly",
                                                 "sub_pc_type": "lu"})

We now isolate the degrees of freedom associated with the vertices and construct a two-level additive Schwarz preconditioner, where the coarse space correction is the inverse of the stiffness matrices associated with the vertices.

def VertexDofs(mesh, fes):
   vertexdofs = BitArray(fes.ndof)
   vertexdofs[:] = False
   for v in mesh.vertices:
      for d in fes.GetDofNrs(v):
         vertexdofs[d] = True
   vertexdofs &= fes.FreeDofs()
   return vertexdofs

vertexdofs = VertexDofs(mesh, fes)
preCoarse = PETScPreconditioner(a.mat, vertexdofs, solverParameters={"pc_type": "hypre"})
pretwo = preCoarse + blockjac
print("-------------------|Additive Schwarz p={}|-------------------".format(order))
gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pretwo, printrates=True)

We can see that the two-level additive Schwarz preconditioner where the coarse space correction is performed using HYPRE is more effective than using just HYPRE for higher-order elements.

Preconditioners performance HYPRE and Two Level Additive Schwarz

Preconditioner

p=1

p=2

p=3

p=4

HYPRE

10 (1.57e-12)

70 (2.22e-13)

42 (3.96e-13)

70 (1.96e-13)

Additive Schwarz

44 (1.96e-12)

45 (1.28e-12)

45 (1.29e-12)

45 (1.45e-12)