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

5.3.2 NGSolve - PETSc interface

We use the ngs2petsc interface to map vectors and matrices between NGSolve and PETSc

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

The Python-module ngsolve.ngs2petsc provides functionality to transfer vectors and matrices between NGSolve and Python.

Make sure that the ipyparallel server can import the module, e.g. by starting the cluster in the current directory.

%%px
import ngsolve.ngs2petsc as n2p
import petsc4py.PETSc as psc

%%px
fes = H1(mesh, order=1, dirichlet="left|bottom")
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

The function CreatePETScMatrix takes an NGSolve matrix, and creates a PETSc matrix from it. A VectorMapping object can map vectors between NGSolve and PETSc.

%%px
psc_mat = n2p.CreatePETScMatrix(a.mat, fes.FreeDofs())
vecmap = n2p.VectorMapping (a.mat.row_pardofs, fes.FreeDofs())

[stderr:0] [d9d8041d895e:15079] Read -1, expected 48512, errno = 1
[stderr:1] [d9d8041d895e:15080] Read -1, expected 33408, errno = 1
[stderr:2] [d9d8041d895e:15081] Read -1, expected 12576, errno = 1

Create PETSc-vectors fitting to the matrix

%%px
psc_f, psc_u = psc_mat.createVecs()
pass

setting up the parallel Krylov-space solver ….

%%px
ksp = psc.KSP()
ksp.create()
ksp.setOperators(psc_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)

moving vectors between NGSolve and PETSc, and solve:

%%px
vecmap.N2P(f.vec, psc_f)
ksp.solve(psc_f, psc_u)
vecmap.P2N(psc_u, gfu.vec)
pass

[stderr:0]
[d9d8041d895e:15079] Read -1, expected 12552, errno = 1
[d9d8041d895e:15079] Read -1, expected 50208, errno = 1
[stderr:1]
[d9d8041d895e:15080] Read -1, expected 9540, errno = 1
[d9d8041d895e:15080] Read -1, expected 16204, errno = 1
[d9d8041d895e:15080] Read -1, expected 64816, errno = 1
[d9d8041d895e:15080] Read -1, expected 38160, errno = 1
[d9d8041d895e:15080] Read -1, expected 5512, errno = 1
[d9d8041d895e:15080] Read -1, expected 6704, errno = 1
[d9d8041d895e:15080] Read -1, expected 6512, errno = 1
[stderr:2]
[d9d8041d895e:15081] Read -1, expected 8204, errno = 1
[d9d8041d895e:15081] Read -1, expected 19168, errno = 1
[d9d8041d895e:15081] Read -1, expected 32816, errno = 1
[d9d8041d895e:15081] Read -1, expected 76672, errno = 1
[d9d8041d895e:15081] Read -1, expected 7144, errno = 1
[d9d8041d895e:15081] Read -1, expected 5720, errno = 1
[d9d8041d895e:15081] Read -1, expected 8832, errno = 1
[stderr:3]
[d9d8041d895e:15082] Read -1, expected 16964, errno = 1
[d9d8041d895e:15082] Read -1, expected 67856, errno = 1
[d9d8041d895e:15082] Read -1, expected 6992, errno = 1

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

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

Draw (gfu[0])

BaseWebGuiScene

PETSc preconditioner for NGSolve

Next we create a PETSc preconditioner, and wrap it into an NGSolve preconditioner:

%%px
a = BilinearForm(grad(u)*grad(v)*dx+u*v*ds)
pre = Preconditioner(a, "gamg")
a.Assemble()
pass

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

[stderr:0]
[d9d8041d895e:15079] Read -1, expected 48512, errno = 1
[d9d8041d895e:15079] Read -1, expected 12552, errno = 1
[d9d8041d895e:15079] Read -1, expected 50208, errno = 1
[stderr:1]
[d9d8041d895e:15080] Read -1, expected 33408, errno = 1
[d9d8041d895e:15080] Read -1, expected 9540, errno = 1
[d9d8041d895e:15080] Read -1, expected 16204, errno = 1
[d9d8041d895e:15080] Read -1, expected 38160, errno = 1
[d9d8041d895e:15080] Read -1, expected 64816, errno = 1
[d9d8041d895e:15080] Read -1, expected 5512, errno = 1
[d9d8041d895e:15080] Read -1, expected 6704, errno = 1
[d9d8041d895e:15080] Read -1, expected 6512, errno = 1
[stderr:2]
[d9d8041d895e:15081] Read -1, expected 12576, errno = 1
[d9d8041d895e:15081] Read -1, expected 8204, errno = 1
[d9d8041d895e:15081] Read -1, expected 19168, errno = 1
[d9d8041d895e:15081] Read -1, expected 76672, errno = 1
[d9d8041d895e:15081] Read -1, expected 32816, errno = 1
[d9d8041d895e:15081] Read -1, expected 7144, errno = 1
[d9d8041d895e:15081] Read -1, expected 5720, errno = 1
[d9d8041d895e:15081] Read -1, expected 8832, errno = 1
[stderr:3]
[d9d8041d895e:15082] Read -1, expected 16964, errno = 1
[d9d8041d895e:15082] Read -1, expected 67856, errno = 1
[d9d8041d895e:15082] Read -1, expected 6992, errno = 1

and use it in an NGSolve - CGSolver:

%%px
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printing=comm.rank==0)
gfu.vec.data = inv * f.vec

[stdout:0]
WARNING: printing is deprecated, use printrates instead!
CG iteration 1, residual = 0.16186373090106865
CG iteration 2, residual = 0.032256334634758486
CG iteration 3, residual = 0.0041873433365436944
CG iteration 4, residual = 0.0005968449121709221
CG iteration 5, residual = 8.461918988055599e-05
CG iteration 6, residual = 1.0980502894984023e-05
CG iteration 7, residual = 1.6480728786485835e-06
CG iteration 8, residual = 2.0479915322333332e-07
CG iteration 9, residual = 2.6152153836992176e-08
CG iteration 10, residual = 3.383295610844336e-09
CG iteration 11, residual = 4.183121875015019e-10
CG iteration 12, residual = 5.703654028140915e-11
CG iteration 13, residual = 7.797932580843257e-12
CG iteration 14, residual = 9.331733543242658e-13
CG iteration 15, residual = 1.2343363523858606e-13

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