PETScUserPreconditioner does not work with PETScKrylovMatrix

Hi,

I came to care about matrix-free in Dolfin, so I decided to have a look at this problem.

The basic problem with the current trunk is obvious from the traceback at the end of the error message:

[0]PETSC ERROR: MatGetFactor() line 3644 in src/mat/interface/matrix.c
[0]PETSC ERROR: PCSetUp_ILU() line 202 in src/ksp/pc/impls/factor/ilu/ilu.c
[0]PETSC ERROR: PCSetUp() line 795 in src/ksp/pc/interface/precon.c

It's in PCSetUp_ILU, but we have told PETScKrylovSolver to use our supplied preconditioner; clearly that's not being registered. Taking a look at dolfin/la/PETScKrylovSolver.cpp, the constructor that takes in a PETScUserPreconditioner assigns it to pc_dolfin, but pc_dolfin is never used again. However, with this patch:

####################
diff -urN dolfin-1.0.0-orig//dolfin/la/PETScKrylovSolver.cpp dolfin-1.0.0//dolfin/la/PETScKrylovSolver.cpp
--- dolfin-1.0.0-orig//dolfin/la/PETScKrylovSolver.cpp 2012-01-28 18:16:39.181840758 +0000
+++ dolfin-1.0.0//dolfin/la/PETScKrylovSolver.cpp 2012-01-28 18:17:13.392732888 +0000
@@ -243,6 +243,11 @@
     preconditioner->set(*this);
     preconditioner_set = true;
   }
+ else if (pc_dolfin && !preconditioner_set)
+ {
+ PETScUserPreconditioner::setup(*_ksp, *pc_dolfin);
+ preconditioner_set = true;
+ }

   // Solve linear system
   if (MPI::process_number() == 0)
####################

the user preconditioner is registered, and it goes into the solve method defined in Python in the example above. However, the solve method as written gives a strange error

    args[0].set_local(args[1])
TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous and uses dtype='d'.

so I replaced it with

####################
from dolfin import *

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Constant(0.0), lambda x, on_boundary: on_boundary)
u = TrialFunction(V); v = TestFunction(V);
A, b = assemble_system( inner(grad(u), grad(v))*dx, Constant(1.0)*v*dx, bc)

class KrylovMatrix(PETScKrylovMatrix) :
    def __init__(self) :
        PETScKrylovMatrix.__init__(self, V.dim(), V.dim())

    def mult(self, *args) :
        y = PETScVector(V.dim())
        A.mult(args[0],y)
        args[1].set_local(y.array())

class IdentityPreconditioner(PETScUserPreconditioner) :
    def __init__(self) :
        PETScUserPreconditioner.__init__(self)

    def solve(self, *args) :
        for i in range(len(args[0])):
          args[0][i] = args[1][i]

y = Function(V)
solve(A,y.vector(),b)
x_petsc = PETScVector(V.dim())

MyPrecon = IdentityPreconditioner()
KrylovSolver = PETScKrylovSolver("cg", MyPrecon)

KrylovSolver.solve(KrylovMatrix(), x_petsc, down_cast(b))
#KrylovSolver.solve(A, x_petsc, down_cast(b)) # works

print (x_petsc - y.vector()).norm("l2")
####################

which now works:

[pef@aislinn:~/src/dolfin/userpc]$ python userpc.py
5.01929900219e-13

Hi. Any word on the above patch?

On 13 February 2012 19:58, Patrick Farrell <email address hidden> wrote:
> Hi. Any word on the above patch?
>

I applied the patch against lp:dolfin (the dev versions), but I can't
run the example code. I get

 Traceback (most recent call last):
  File "test.py", line 33, in <module>
    KrylovSolver.solve(KrylovMatrix(), x_petsc, down_cast(b))
  File "test.py", line 23, in solve
    for i in range(len(args[0])):
TypeError: object of type 'SwigPyObject' has no len()

Could you try getting the demo code to run against lp;dolfin?

Garth

> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/768233
>
> Title:
>  PETScUserPreconditioner does not work with PETScKrylovMatrix
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/768233/+subscriptions

Thanks for the patch and sorry for the long time taken to get it applied. I am working on a unit test for this. Will push when all looks fine.

I have now applied the patch, and added a unittest for MatrixFree krylov
solver. The interface of PETScKrylovSolver have to be improved as it is
no problem to break the solve process by for example not providing a
Userdefined Preconditioner.

Johan

Great, thanks!