NGSolve - PETSc interface

107. NGSolve - PETSc interface#

We use the ngs2petsc interface (shipped with NGSolve) to map vectors and matrices between NGSolve and PETSc

An extended interface to PETSc is available by Stefano Zampini and Umberto Zerbinati from https://ngspetsc.readthedocs.io/en/latest/ngsPETSc.html

from ipyparallel import Cluster
c = await Cluster(engines="mpi").start_and_connect(n=4, activate=True)
c.ids

Starting 4 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>

[0, 1, 2, 3]

%%px
from ngsolve import *
comm = MPI.COMM_WORLD
ngmesh = unit_square.GenerateMesh(maxh=0.1, comm=comm)
    
for l in range(3):
    ngmesh.Refine()
mesh = Mesh(ngmesh)

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

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

%%px
# psc_mat.view()

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

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

Draw (gfu[0]);

107.1. 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();

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.15523799534158494     
CG iteration 2, residual = 0.030515901578059508     
CG iteration 3, residual = 0.005254427671245441     
CG iteration 4, residual = 0.000991222354408726     
CG iteration 5, residual = 0.00017826167799212423     
CG iteration 6, residual = 3.6929115679800176e-05     
CG iteration 7, residual = 6.557805237886243e-06     
CG iteration 8, residual = 1.2378725412901232e-06     
CG iteration 9, residual = 2.2086097459486483e-07     
CG iteration 10, residual = 3.782341314152293e-08     
CG iteration 11, residual = 6.357967349566024e-09     
CG iteration 12, residual = 1.0292339438117692e-09     
CG iteration 13, residual = 1.7201887732936172e-10     
CG iteration 14, residual = 2.643158606468948e-11     
CG iteration 15, residual = 4.314906419146038e-12     
CG iteration 16, residual = 6.951824746598483e-13     
CG iteration 17, residual = 1.1093622681085737e-13     

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

c.shutdown(hub=True)
Contents