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 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.occ import unit_square
comm = MPI.COMM_WORLD

ngmesh = unit_square.GenerateMesh(maxh=0.1, comm=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

[stderr:2] [e278b3a455ef:67171] Read -1, expected 4096, errno = 1

[stderr:1] [e278b3a455ef:67170] Read -1, expected 4096, errno = 1

[stderr:3] [e278b3a455ef:67172] Read -1, expected 4096, errno = 1

%%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] [e278b3a455ef:67171] Read -1, expected 15168, errno = 1

[stderr:0] [e278b3a455ef:67169] Read -1, expected 48832, errno = 1

[stderr:1] [e278b3a455ef:67170] Read -1, expected 35488, errno = 1

Create PETSc-vectors fitting to the matrix

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

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

[stderr:2] [e278b3a455ef:67171] Read -1, expected 9904, errno = 1
[e278b3a455ef:67171] Read -1, expected 18176, errno = 1
[e278b3a455ef:67171] Read -1, expected 39616, errno = 1
[e278b3a455ef:67171] Read -1, expected 72704, errno = 1
[e278b3a455ef:67171] Read -1, expected 8072, errno = 1
[e278b3a455ef:67171] Read -1, expected 5904, errno = 1
[e278b3a455ef:67171] Read -1, expected 8416, errno = 1

[stderr:3] [e278b3a455ef:67172] Read -1, expected 19700, errno = 1
[e278b3a455ef:67172] Read -1, expected 78800, errno = 1
[e278b3a455ef:67172] Read -1, expected 4352, errno = 1
[e278b3a455ef:67172] Read -1, expected 7392, errno = 1

[stderr:1] [e278b3a455ef:67170] Read -1, expected 8672, errno = 1
[e278b3a455ef:67170] Read -1, expected 16948, errno = 1
[e278b3a455ef:67170] Read -1, expected 34688, errno = 1
[e278b3a455ef:67170] Read -1, expected 67792, errno = 1
[e278b3a455ef:67170] Read -1, expected 4432, errno = 1
[e278b3a455ef:67170] Read -1, expected 5040, errno = 1
[e278b3a455ef:67170] Read -1, expected 6104, errno = 1
[e278b3a455ef:67170] Read -1, expected 5776, errno = 1

[stderr:0] [e278b3a455ef:67169] Read -1, expected 11384, errno = 1
[e278b3a455ef:67169] Read -1, expected 45536, 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, "petsc", pctype="gamg", levels=10)
pre = Preconditioner(a, "gamg")
a.Assemble();

[stderr:1] [e278b3a455ef:67170] Read -1, expected 35488, errno = 1
[e278b3a455ef:67170] Read -1, expected 8672, errno = 1
[e278b3a455ef:67170] Read -1, expected 16948, errno = 1
[e278b3a455ef:67170] Read -1, expected 34688, errno = 1
[e278b3a455ef:67170] Read -1, expected 67792, errno = 1
[e278b3a455ef:67170] Read -1, expected 4432, errno = 1
[e278b3a455ef:67170] Read -1, expected 5040, errno = 1
[e278b3a455ef:67170] Read -1, expected 6104, errno = 1
[e278b3a455ef:67170] Read -1, expected 5776, errno = 1

[stderr:3] [e278b3a455ef:67172] Read -1, expected 19700, errno = 1
[e278b3a455ef:67172] Read -1, expected 78800, errno = 1
[e278b3a455ef:67172] Read -1, expected 4352, errno = 1
[e278b3a455ef:67172] Read -1, expected 7392, errno = 1

[stderr:2] [e278b3a455ef:67171] Read -1, expected 15168, errno = 1
[e278b3a455ef:67171] Read -1, expected 9904, errno = 1
[e278b3a455ef:67171] Read -1, expected 18176, errno = 1
[e278b3a455ef:67171] Read -1, expected 39616, errno = 1
[e278b3a455ef:67171] Read -1, expected 72704, errno = 1
[e278b3a455ef:67171] Read -1, expected 8072, errno = 1
[e278b3a455ef:67171] Read -1, expected 5904, errno = 1
[e278b3a455ef:67171] Read -1, expected 8416, errno = 1

[stderr:0] [e278b3a455ef:67169] Read -1, expected 48832, errno = 1
[e278b3a455ef:67169] Read -1, expected 11384, errno = 1
[e278b3a455ef:67169] Read -1, expected 45536, errno = 1

and use it in an NGSolve - CGSolver:

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

[stdout:0] CG iteration 1, residual = 0.1616451275202024
CG iteration 2, residual = 0.03171900063988554
CG iteration 3, residual = 0.004198128230263991
CG iteration 4, residual = 0.0006197236145500212
CG iteration 5, residual = 8.613269419324584e-05
CG iteration 6, residual = 1.0597037303146737e-05
CG iteration 7, residual = 1.5632460635603687e-06
CG iteration 8, residual = 1.85246537964527e-07
CG iteration 9, residual = 2.4048716218910356e-08
CG iteration 10, residual = 3.6369768538179375e-09
CG iteration 11, residual = 5.028733998453997e-10
CG iteration 12, residual = 7.216444708125522e-11
CG iteration 13, residual = 9.633723762360763e-12
CG iteration 14, residual = 1.133884901416395e-12
CG iteration 15, residual = 1.5325405389402187e-13

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