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)

[stdout:0]
 Generate Mesh from spline geometry
 Boundary mesh done, np = 40
 CalcLocalH: 40 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 call metis 5 ...
 metis start
 metis complete
 Send/Receive mesh
 Update mesh topology
 Sending nr of elements
 Building vertex/proc mapping
 Sending Vertices - vertices
 Sending Vertices - identifications
 Sending Vertices - distprocs
 Sending elements
 Sending Face Descriptors
 Sending Surface elements
 Sending Edge Segments
 Point-Elements ...
 now wait ...
 Sending names
 wait for names
 Clean up local memory
 Update mesh topology
 send mesh complete
 update parallel topology
 Refine mesh
 update parallel topology
 Update mesh topology
 Update clusters
 update parallel topology
 Refine mesh
 update parallel topology
 Update mesh topology
 Update clusters
 update parallel topology
 Update clusters
 update parallel topology
[stdout:1]  p1: got 0 elements and 76 surface elements
[stdout:2]  p2: got 0 elements and 81 surface elements
[stdout:3]  p3: got 0 elements and 75 surface elements

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

[stdout:0]
assemble VOL element 0/0
assemble BND element 0/0
assemble VOL element 0/0

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:0] [d9d8041d895e:15079] Read -1, expected 50640, errno = 1
[stderr:1] [d9d8041d895e:15080] Read -1, expected 32064, errno = 1
[stderr:2] [d9d8041d895e:15081] Read -1, expected 13008, errno = 1

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

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

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

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

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

[stderr:0]
[d9d8041d895e:15079] Read -1, expected 15344, errno = 1
[d9d8041d895e:15079] Read -1, expected 61376, errno = 1
[d9d8041d895e:15079] Read -1, expected 4736, errno = 1
[stderr:1]
[d9d8041d895e:15080] Read -1, expected 11536, errno = 1
[d9d8041d895e:15080] Read -1, expected 15788, errno = 1
[d9d8041d895e:15080] Read -1, expected 46144, errno = 1
[d9d8041d895e:15080] Read -1, expected 63152, errno = 1
[d9d8041d895e:15080] Read -1, expected 4096, errno = 1
[d9d8041d895e:15080] Read -1, expected 6112, errno = 1
[d9d8041d895e:15080] Read -1, expected 6968, errno = 1
[d9d8041d895e:15080] Read -1, expected 6880, errno = 1
[stderr:2]
[d9d8041d895e:15081] Read -1, expected 8580, errno = 1
[d9d8041d895e:15081] Read -1, expected 19656, errno = 1
[d9d8041d895e:15081] Read -1, expected 34320, errno = 1
[d9d8041d895e:15081] Read -1, expected 78624, errno = 1
[d9d8041d895e:15081] Read -1, expected 7400, errno = 1
[d9d8041d895e:15081] Read -1, expected 5392, errno = 1
[d9d8041d895e:15081] Read -1, expected 9360, errno = 1
[stderr:3]
[d9d8041d895e:15082] Read -1, expected 16972, errno = 1
[d9d8041d895e:15082] Read -1, expected 67888, errno = 1
[d9d8041d895e:15082] Read -1, expected 7696, 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])

 Update mesh topology
 Update clusters
 Update mesh topology
 Update clusters
 Update mesh topology
 Update clusters
 Update mesh topology
 Update clusters

BaseWebGuiScene