petsc linking

hello, I would like to know if its possible to use PETSc solvers in NGsolve.
Regards
Tom

Hi Tom,
actually I was interested in this for some time now, your question pushed me to try it
It should be possible using the petsc4py interface. Note that I have no experience at all with petsc and the documentation of their Python interface is not existent… Maybe some things can be done more elegant
I’ve modified the poisson.py example to use a petsc solver, interfacing with numpy/scipy. Note that on our side no copies are made only type conversions - petsc may do copy / mpi communcations

First we set up the ngsolve problem:

``````
import petsc4py
petsc4py.init()
import petsc4py.PETSc as psc
import numpy as np
from ngsolve import *
from scipy.sparse import csr_matrix
from netgen.geom2d import unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = H1(mesh, order=3)
u,v = fes.TnT()

f = LinearForm(fes)
f += SymbolicLFI(32 * (y*(1-y)+x*(1-x)) * v)

a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v))
a += SymbolicBFI(1e6* u*v, definedon=mesh.Boundaries("top|bottom|right|left"))

a.Assemble()
f.Assemble()

gfu = GridFunction(fes)``````

Then wrap our data with petsc classes and solve there (and draw the solution):

``````fpsc = psc.Vec()
fpsc.createWithArray(f.vec.FV().NumPy())

spy_csr = csr_matrix(a.mat.CSR(), shape=(a.mat.height, a.mat.width))

apsc = psc.Mat().createAIJ(size=spy_csr.shape)
apsc.setUp()
(Istart, Iend) = apsc.getOwnershipRange()
ai = spy_csr.indptr[Istart:Iend+1] - spy_csr.indptr[Istart]
aj = spy_csr.indices[spy_csr.indptr[Istart]:spy_csr.indptr[Iend]]
av = spy_csr.data[spy_csr.indptr[Istart]:spy_csr.indptr[Iend]]

apsc.setValuesCSR(ai,aj,av)
apsc.assemble()

ksp = psc.KSP()
ksp.create()
ksp.setOperators(apsc)

upsc = psc.Vec()
upsc.createWithArray(gfu.vec.FV().NumPy())

ksp.solve(fpsc,upsc)

Draw (gfu)``````

Since petsc writes the solution into our vector I assume at least for the vectors there are no copies involved.

I don’t know how to set dirichlet (non free dofs) for petsc so that it only solves on sub matrices, thats why I implemented the dirichlet boundaries with penalty terms.

Best
Christopher

Attachment: ngsolve_petsc.py

Thank you very much. I’m would like to try my problem in NGSolve as it supports p-refinement. I work on Maxwell’s equations and will get back to you after trying my problem using this script. Thanks again and merry x’mas !
regards