104. Using PETSc

PETSc https://www.mcs.anl.gov/petsc
Portable, Extensible Toolkit for Scientific Computation Toolkit

PETSc is the standard library for parallel solvers. Many other libraries offer interfaces to PETSc.

We learn how to interface PETSc from NGSolve.

from ipyparallel import Cluster
c = await Cluster(engines="mpi").start_and_connect(n=4, activate=True)
c.ids

Starting 4 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>

[0, 1, 2, 3]

%%px
from ngsolve import *
comm = MPI.COMM_WORLD
ngmesh = unit_square.GenerateMesh(maxh=0.1, comm=comm)
    
for l in range(3):
    ngmesh.Refine()
mesh = Mesh(ngmesh)

We import PETSc via its Python interface petsc4py:

%%px
import numpy as np
import petsc4py.PETSc as psc

%%px
fes = H1(mesh, order=1)
print ("ndof =", fes.ndof, "ndglob =", fes.ndofglobal)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx+u*v*ds).Assemble()
f = LinearForm(x*v*dx).Assemble()
gfu = GridFunction(fes)

[stdout:0] ndof = 0 ndglob = 7585

[stdout:3] ndof = 2481 ndglob = 7585

[stdout:1] ndof = 2593 ndglob = 7585

[stdout:2] ndof = 2657 ndglob = 7585

We convert the local sparse matrix to a local PETSc AIJ matrix:

%%px
locmat = a.mat.local_mat
val,col,ind = locmat.CSR()
ind = np.array(ind, dtype='int32')
apsc_loc = psc.Mat().createAIJ(size=(locmat.height, locmat.width), csr=(ind,col,val), comm=MPI.COMM_SELF)

The NGSolve ParallelDofs object corresponds to a PETSc IndexSet object, which we create next. In PETSc, dofs are globally enumerated, what is not the case in NGSolve. For this purpose, a ParallelDofs class can generate a globally consistent enumeration of dofs. The generated globnums array contains the global dof-numbers of the local dofs.

%%px
pardofs = fes.ParallelDofs()
globnums, nglob = pardofs.EnumerateGlobally() 
# print (list(globnums))

iset = psc.IS().createGeneral (indices=globnums, comm=comm)
lgmap = psc.LGMap().createIS(iset)

We can now create the global matrix using our local2global map:

%%px
mat = psc.Mat().create(comm=comm)
mat.setSizes(size=(nglob,nglob), bsize=1)
mat.setType(psc.Mat.Type.IS)
mat.setLGMap(lgmap)
mat.setISLocalMat(apsc_loc)
mat.assemble()
mat.convert("mpiaij")

Out[2:6]: <petsc4py.PETSc.Mat at 0x10a18d4e0>

Out[3:6]: <petsc4py.PETSc.Mat at 0x11c921620>

Out[0:6]: <petsc4py.PETSc.Mat at 0x1166e1670>

Out[1:6]: <petsc4py.PETSc.Mat at 0x119738ef0>

We copy the right-hand side vector to a PETSc vector. PETSc vectors (knowing about the global connectivity) are created from the parallel PETSc matrix. We copy the NGSolve-vector to a local sub-vector:

%%px

f.vec.Cumulate()

v1, v2 = mat.createVecs()

v2loc = v2.getSubVector(iset)
v2loc.getArray()[:] = f.vec 
v2.restoreSubVector(iset, v2loc)

%%px
ksp = psc.KSP()
ksp.create()
ksp.setOperators(mat)
ksp.setType(psc.KSP.Type.CG)
ksp.setNormType(psc.KSP.NormType.NORM_NATURAL)
ksp.getPC().setType("gamg")
ksp.setTolerances(rtol=1e-6, atol=0, divtol=1e16, max_it=400)
ksp.solve(v2,v1)   

print ("petsc-its =", ksp.its)

[stdout:0] petsc-its = 8

[stdout:2] petsc-its = 8

[stdout:1] petsc-its = 8

[stdout:3] petsc-its = 8

%%px
v1loc = v1.getSubVector(iset)
gfu.vec[:] = v1loc.getArray()

%%px
print ("<f,u> =", InnerProduct(f.vec, gfu.vec))

[stdout:2] <f,u> = 0.07795637134606541

[stdout:0] <f,u> = 0.07795637134606541

[stdout:1] <f,u> = 0.07795637134606541

[stdout:3] <f,u> = 0.07795637134606541

from ngsolve.webgui import Draw
gfu = c[:]["gfu"]

Draw (gfu[0], order=1);

c.shutdown(hub=True)