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 Client
c = Client(profile='mpi')
c.ids

[0, 1, 2, 3]

%%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")

Out[0:26]: <petsc4py.PETSc.Mat at 0x7f74caedd770>

Out[1:26]: <petsc4py.PETSc.Mat at 0x7f6aa7c31ae0>

Out[2:26]: <petsc4py.PETSc.Mat at 0x7f0e17fd9040>

Out[3:26]: <petsc4py.PETSc.Mat at 0x7f5511c38040>

%%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:0] petsc-its = 7
[stdout:1] petsc-its = 7
[stdout:2] petsc-its = 7
[stdout:3] petsc-its = 7

%%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