- Thank you received: 0
petsc linking
5 years 11 months ago - 5 years 11 months ago #1349
by christopher
Replied by christopher on topic petsc linking
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:
Then wrap our data with petsc classes and solve there (and draw the solution):
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
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:
Code:
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):
Code:
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
Attachments:
Time to create page: 0.098 seconds