# 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)