petsc linking

Main Forum Old-Forum - do not post here
1

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

2

Hi Tom,
actually I was interested in this for some time now, your question pushed me to try it :wink:
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 :wink:
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

3

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