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)

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)

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:2] [28b404a97ea5:13863] Read -1, expected 18256, errno = 1

[stderr:0] [28b404a97ea5:13861] Read -1, expected 48352, errno = 1

[stderr:1] [28b404a97ea5:13862] Read -1, expected 36240, 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:1] [28b404a97ea5:13862] Read -1, expected 7992, errno = 1
[28b404a97ea5:13862] Read -1, expected 16864, errno = 1
[28b404a97ea5:13862] Read -1, expected 31968, errno = 1
[28b404a97ea5:13862] Read -1, expected 67456, errno = 1
[28b404a97ea5:13862] Read -1, expected 4800, errno = 1
[28b404a97ea5:13862] Read -1, expected 4424, errno = 1
[28b404a97ea5:13862] Read -1, expected 5288, errno = 1
[28b404a97ea5:13862] Read -1, expected 4752, errno = 1

[stderr:0] [28b404a97ea5:13861] Read -1, expected 10588, errno = 1
[28b404a97ea5:13861] Read -1, expected 42352, errno = 1

[stderr:2] [28b404a97ea5:13863] Read -1, expected 11692, errno = 1
[28b404a97ea5:13863] Read -1, expected 16364, errno = 1
[28b404a97ea5:13863] Read -1, expected 46768, errno = 1
[28b404a97ea5:13863] Read -1, expected 65456, errno = 1
[28b404a97ea5:13863] Read -1, expected 4480, errno = 1
[28b404a97ea5:13863] Read -1, expected 8960, errno = 1
[28b404a97ea5:13863] Read -1, expected 6056, errno = 1
[28b404a97ea5:13863] Read -1, expected 8464, errno = 1

[stderr:3] [28b404a97ea5:13864] Read -1, expected 22384, errno = 1
[28b404a97ea5:13864] Read -1, expected 89536, errno = 1
[28b404a97ea5:13864] Read -1, expected 8848, errno = 1

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

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

[stderr:0] [28b404a97ea5:13861] Read -1, expected 48352, errno = 1
[28b404a97ea5:13861] Read -1, expected 10588, errno = 1
[28b404a97ea5:13861] Read -1, expected 42352, errno = 1

[stderr:2] [28b404a97ea5:13863] Read -1, expected 18256, errno = 1
[28b404a97ea5:13863] Read -1, expected 11692, errno = 1
[28b404a97ea5:13863] Read -1, expected 16364, errno = 1
[28b404a97ea5:13863] Read -1, expected 65456, errno = 1
[28b404a97ea5:13863] Read -1, expected 46768, errno = 1
[28b404a97ea5:13863] Read -1, expected 4480, errno = 1
[28b404a97ea5:13863] Read -1, expected 8960, errno = 1
[28b404a97ea5:13863] Read -1, expected 6056, errno = 1
[28b404a97ea5:13863] Read -1, expected 8464, errno = 1

[stderr:3] [28b404a97ea5:13864] Read -1, expected 22384, errno = 1
[28b404a97ea5:13864] Read -1, expected 89536, errno = 1
[28b404a97ea5:13864] Read -1, expected 8848, errno = 1

[stderr:1] [28b404a97ea5:13862] Read -1, expected 36240, errno = 1
[28b404a97ea5:13862] Read -1, expected 7992, errno = 1
[28b404a97ea5:13862] Read -1, expected 16864, errno = 1
[28b404a97ea5:13862] Read -1, expected 31968, errno = 1
[28b404a97ea5:13862] Read -1, expected 67456, errno = 1
[28b404a97ea5:13862] Read -1, expected 4800, errno = 1
[28b404a97ea5:13862] Read -1, expected 4424, errno = 1
[28b404a97ea5:13862] Read -1, expected 5288, errno = 1
[28b404a97ea5:13862] Read -1, expected 4752, 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.1618787001637246
CG iteration 2, residual = 0.031964365197057684
CG iteration 3, residual = 0.0038382478132801393
CG iteration 4, residual = 0.0005682457460902522
CG iteration 5, residual = 7.955315146595451e-05
CG iteration 6, residual = 9.699739268709486e-06
CG iteration 7, residual = 1.5069797657401885e-06
CG iteration 8, residual = 1.9165943622363737e-07
CG iteration 9, residual = 2.1994160499747156e-08
CG iteration 10, residual = 3.0749779591135743e-09
CG iteration 11, residual = 4.036559619891373e-10
CG iteration 12, residual = 5.060126332469917e-11
CG iteration 13, residual = 6.5441849062792505e-12
CG iteration 14, residual = 8.105131770777225e-13
CG iteration 15, residual = 1.054169477629426e-13

gfu = c[:]["gfu"]
Draw (gfu[0])

BaseWebGuiScene