Matrix-free approach

Linear Algebra
1

Hi everyone,

I would like to solve a Poisson problem with a matrix-free approach.
For now, I managed to do it by using a PETSc ksp solver and defining a mult routine for which I used the ufl.action function (see example below).

However, I have the following problems:

  1. This example does not work in parallel.
  2. The number of iterations of the conjugate gradient is too important. How can I define a matrix-free preconditioner for this problem ?

Cheers,
Fabien

from dolfin import *
import petsc4py.PETSc as petsc

n = 100
mesh = UnitSquareMesh(n,n)

class Laplace:
    def __init__(self, mesh):
        self.V = FunctionSpace(mesh,"CG",1)
        u = TrialFunction(self.V)
        v = TestFunction(self.V)
        self.a = inner(grad(u),grad(v))*dx
        self.L = inner(Constant(1.),v)*dx
        def bounds(x, on_boundary):
            return on_boundary
        self.bc = DirichletBC(self.V,Constant(0.),bounds)
        self.sol = Function(self.V)   
        self.b = PETScVector()
        self.A = PETScMatrix()
        assemble(self.L, tensor=self.b)
        self.bc.apply(self.b)

    def mult(self, mat, X, Y):
        csol = self.sol.vector().vec().getArray(readonly=0)
        XX = X.getArray(readonly=1)
        csol[:] = XX[:]
        Act = action(self.a, self.sol)
        AX = PETScVector()
        assemble(Act, tensor=AX)
        self.bc.apply(AX)
        Y.getArray(readonly=0)[:] = AX.vec().getArray(readonly=1)[:]
    
    def setup_ksp(self, solver='cg', prec='none'):
        fM = assemble(TrialFunction(self.V)*TestFunction(self.V)*dx)
        M = petsc.Mat()
        M.create(comm=petsc.COMM_WORLD)
        M.setSizes(as_backend_type(fM).mat().getSizes())
        M.setType(petsc.Mat.Type.PYTHON)
        M.setPythonContext(self)
        M.setUp()

        PETScOptions.set("ksp_monitor_true_residual")
        self.ksp = petsc.KSP().create()
        self.ksp.setType(solver)
        pc = self.ksp.getPC()
        pc.setType(prec)
        self.ksp.setOperators(M)
        self.ksp.setName("laplace")
        self.ksp.setTolerances(rtol=1e-6)
        self.ksp.setFromOptions()

        x, b = M.getVecs()
        b.getArray(readonly=0)[:] = self.b.vec().getArray(readonly=1)[:]
        return x, b

    def solve(self, B, X):
        self.ksp.solve(B, X)
        self.sol.vector().vec().getArray(readonly=0)[:] = X.getArray(readonly=1)[:]


l = Laplace(mesh)
x, b = l.setup_ksp()
l.solve(b,x)
comm = MPI.comm_world
with XDMFFile(comm,"laplace_matrixfree.xdmf") as infile:
    infile.write(l.sol)