Comment 1 for bug 768233

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