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

%%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:0] [5c736296b3b3:69860] Read -1, expected 48864, errno = 1

[stderr:1] [5c736296b3b3:69861] Read -1, expected 33824, errno = 1

[stderr:2] [5c736296b3b3:69862] Read -1, expected 17328, 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:1] [5c736296b3b3:69861] Read -1, expected 9512, errno = 1
[5c736296b3b3:69861] Read -1, expected 15980, errno = 1
[5c736296b3b3:69861] Read -1, expected 38048, errno = 1
[5c736296b3b3:69861] Read -1, expected 63920, errno = 1
[5c736296b3b3:69861] Read -1, expected 5200, errno = 1
[5c736296b3b3:69861] Read -1, expected 5672, errno = 1
[5c736296b3b3:69861] Read -1, expected 5680, errno = 1

[stderr:3] [5c736296b3b3:69863] Read -1, expected 22316, errno = 1
[5c736296b3b3:69863] Read -1, expected 89264, errno = 1
[5c736296b3b3:69863] Read -1, expected 8032, errno = 1

[stderr:0] [5c736296b3b3:69860] Read -1, expected 12480, errno = 1
[5c736296b3b3:69860] Read -1, expected 49920, errno = 1

[stderr:2] [5c736296b3b3:69862] Read -1, expected 11312, errno = 1
[5c736296b3b3:69862] Read -1, expected 16740, errno = 1
[5c736296b3b3:69862] Read -1, expected 45248, errno = 1
[5c736296b3b3:69862] Read -1, expected 66960, errno = 1
[5c736296b3b3:69862] Read -1, expected 4392, errno = 1
[5c736296b3b3:69862] Read -1, expected 8784, errno = 1
[5c736296b3b3:69862] Read -1, expected 5640, errno = 1
[5c736296b3b3:69862] Read -1, expected 7664, 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:0] [5c736296b3b3:69860] Read -1, expected 48864, errno = 1
[5c736296b3b3:69860] Read -1, expected 12480, errno = 1
[5c736296b3b3:69860] Read -1, expected 49920, errno = 1

[stderr:1] [5c736296b3b3:69861] Read -1, expected 33824, errno = 1
[5c736296b3b3:69861] Read -1, expected 9512, errno = 1
[5c736296b3b3:69861] Read -1, expected 15980, errno = 1
[5c736296b3b3:69861] Read -1, expected 38048, errno = 1
[5c736296b3b3:69861] Read -1, expected 63920, errno = 1
[5c736296b3b3:69861] Read -1, expected 5200, errno = 1
[5c736296b3b3:69861] Read -1, expected 5672, errno = 1
[5c736296b3b3:69861] Read -1, expected 5680, errno = 1

[stderr:3] [5c736296b3b3:69863] Read -1, expected 22316, errno = 1
[5c736296b3b3:69863] Read -1, expected 89264, errno = 1
[5c736296b3b3:69863] Read -1, expected 8032, errno = 1

[stderr:2] [5c736296b3b3:69862] Read -1, expected 17328, errno = 1
[5c736296b3b3:69862] Read -1, expected 11312, errno = 1
[5c736296b3b3:69862] Read -1, expected 16740, errno = 1
[5c736296b3b3:69862] Read -1, expected 66960, errno = 1
[5c736296b3b3:69862] Read -1, expected 45248, errno = 1
[5c736296b3b3:69862] Read -1, expected 4392, errno = 1
[5c736296b3b3:69862] Read -1, expected 8784, errno = 1
[5c736296b3b3:69862] Read -1, expected 5640, errno = 1
[5c736296b3b3:69862] Read -1, expected 7664, 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.1612701068617255
CG iteration 2, residual = 0.031430324727711303
CG iteration 3, residual = 0.004236806109474581
CG iteration 4, residual = 0.0006735602561428321
CG iteration 5, residual = 9.254951509649878e-05
CG iteration 6, residual = 1.1543042891878422e-05
CG iteration 7, residual = 1.6313785246666158e-06
CG iteration 8, residual = 2.0609628931866596e-07
CG iteration 9, residual = 2.8655513283897342e-08
CG iteration 10, residual = 4.197371013467841e-09
CG iteration 11, residual = 4.948642085120869e-10
CG iteration 12, residual = 6.073045364589212e-11
CG iteration 13, residual = 8.76125554964775e-12
CG iteration 14, residual = 1.1080631830811685e-12
CG iteration 15, residual = 1.4458991596493695e-13

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