This page was generated from unit-5a.3-petsc/petsc.ipynb.

5.3.1 Using PETSc

We learn how we can interface the popular parallel toolkit PETSc for solving linear equations and much more. We use the Python interface petsc4py. The whole Python file is n2p_ex1.py

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

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

%%px
from ngsolve import *
from netgen.geom2d import unit_square
comm = MPI.COMM_WORLD
if comm.rank == 0:
    ngmesh = unit_square.GenerateMesh(maxh=0.1).Distribute(comm)
else:
    ngmesh = netgen.meshing.Mesh.Receive(comm)

for l in range(2):
    ngmesh.Refine()
mesh = Mesh(ngmesh)

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

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

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().createPython(size=nglob, comm=comm)
mat.setType(psc.Mat.Type.IS)
mat.setLGMap(lgmap)
mat.setISLocalMat(apsc_loc)
mat.assemble()
mat.convert("mpiaij")

[stderr:1] [d8c32f2f330d:70985] Read -1, expected 34480, errno = 1

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

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

[stderr:2] [d8c32f2f330d:70986] Read -1, expected 16400, errno = 1

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

[stderr:0] [d8c32f2f330d:70984] Read -1, expected 50576, errno = 1

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

%%px

f.vec.Cumulate()

v1, v2 = mat.createVecs()

v2loc = v2.getSubVector(iset)
v2loc.getArray()[:] = f.vec.FV()
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:1] petsc-its = 7

[stdout:3] petsc-its = 7

[stdout:0] petsc-its = 7

[stdout:2] petsc-its = 7

[stderr:1] [d8c32f2f330d:70985] Read -1, expected 9960, errno = 1
[d8c32f2f330d:70985] Read -1, expected 16508, errno = 1
[d8c32f2f330d:70985] Read -1, expected 39840, errno = 1
[d8c32f2f330d:70985] Read -1, expected 66032, errno = 1
[d8c32f2f330d:70985] Read -1, expected 4320, errno = 1
[d8c32f2f330d:70985] Read -1, expected 5504, errno = 1
[d8c32f2f330d:70985] Read -1, expected 6216, errno = 1
[d8c32f2f330d:70985] Read -1, expected 6336, errno = 1

[stderr:3] [d8c32f2f330d:70987] Read -1, expected 21272, errno = 1
[d8c32f2f330d:70987] Read -1, expected 85088, errno = 1
[d8c32f2f330d:70987] Read -1, expected 8368, errno = 1

[stderr:0] [d8c32f2f330d:70984] Read -1, expected 13156, errno = 1
[d8c32f2f330d:70984] Read -1, expected 52624, errno = 1

[stderr:2] [d8c32f2f330d:70986] Read -1, expected 10804, errno = 1
[d8c32f2f330d:70986] Read -1, expected 18316, errno = 1
[d8c32f2f330d:70986] Read -1, expected 73264, errno = 1
[d8c32f2f330d:70986] Read -1, expected 43216, errno = 1
[d8c32f2f330d:70986] Read -1, expected 4336, errno = 1
[d8c32f2f330d:70986] Read -1, expected 8672, errno = 1
[d8c32f2f330d:70986] Read -1, expected 5816, errno = 1
[d8c32f2f330d:70986] Read -1, expected 8368, errno = 1

%%px
v1loc = v1.getSubVector(iset)
for i in range(len(gfu.vec)):
    gfu.vec.FV()[i] = v1loc.getArray()[i]

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

BaseWebGuiScene