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

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

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

Draw (gfu[0])

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

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]
iteration 0 error = 0.16186373090106865
iteration 1 error = 0.03225633463475834
iteration 2 error = 0.00418734333654368
iteration 3 error = 0.0005968449121709168
iteration 4 error = 8.46191898805545e-05
iteration 5 error = 1.0980502894984626e-05
iteration 6 error = 1.6480728786493365e-06
iteration 7 error = 2.0479915322434306e-07
iteration 8 error = 2.615215383748639e-08
iteration 9 error = 3.3832956106055766e-09
iteration 10 error = 4.183121870187323e-10
iteration 11 error = 5.703654006011637e-11
iteration 12 error = 7.79793339014583e-12
iteration 13 error = 9.331741996768944e-13
iteration 14 error = 1.2343443410749943e-13

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