PETScKrylovSolver won't solve in version 1.4.0

Question

The following code seems to work in 1.6.0 but gives an error in 1.4.0
(see also https://github.com/Juanlu001/fenics-recipes/issues/28)
Any pointers to where something could be wrong?

from dolfin import *
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))* dx
f = Expression('x[0]*x[1]')
b = f * v * dx
p = a + u * v * dx

A = assemble(a)
B = assemble(b)
P = assemble(p)

class CompositeOperator(object):
    def __init__(self, A, B, C):
        self.A_mat = as_backend_type(A).mat()
        self.B_vec = as_backend_type(B).vec()
        self.C_vec = as_backend_type(C).vec()

    def mult(self, mat, x, y):
        self.A_mat.mult(x, y)
        a = self.C_vec.dot(x)
        y.axpy(a, self.B_vec)

    def as_petscmatrix(self):
        from petsc4py import PETSc
        mat = PETSc.Mat().createPython(self.A_mat.getSizes(), comm = mpi_comm_world())
        mat.setPythonContext(self)
        return PETScMatrix(mat)

# Apply boundary conditions
bc = DirichletBC(V, 1., "x[1]<DOLFIN_EPS")
bc.apply(A)
bc.apply(P)
C = B.copy()
bc.homogenize()
bc.apply(B)


A_B = CompositeOperator(A, B, C).as_petscmatrix()

solver = PETScKrylovSolver()
solver.set_operators(A_B, P)

solution = Function(V)
solver.solve(solution.vector(), C)

The error:

Traceback (most recent call last):
  File "<path>/temp.py", line 50, in <module>
    solver.solve(solution.vector(), C)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 73.
*** Where:   This error was encountered inside /home/juanlu/miniconda/conda-bld/work/dolfin-1.4.0/dolfin/la/PETScKrylovSolver.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  
*** -------------------------------------------------------------------------

Error 73 is "object in argument is in wrong state, e.g. unassembled mat".
http://www.mcs.anl.gov/petsc/petsc-master/include/petscerror.h.html

Øyvind Evju · Answer 1 · Sep 9, 2015

Best answer

Try finalizing assembly of the matrix by doing

A_B.apply("add")

or

A_B.mat().assemble()

before solve.

(I would advice you to use the most recent version of dolfin.)

answered Sep 9, 2015 by Øyvind Evju FEniCS Expert (17,700 points)
selected Sep 11, 2015 by maartent

Both give an error, eg. for A_B.mat().assemble()

Traceback (most recent call last):
  File "<path>/temp.py", line 42, in <module>
    A_B.mat().assemble()
  File "PETSc/Mat.pyx", line 911, in petsc4py.PETSc.Mat.assemble (src/petsc4py.PETSc.c:104649)
petsc4py.PETSc.Error: error code 73
[0] MatAssemblyBegin() line 4794 in src/mat/interface/matrix.c
[0] Object is in wrong state
[0] Must call MatXXXSetPreallocation() or MatSetUp() on argument 1 "mat" before MatAssemblyBegin()

What is meant by 'assembling' the matrix in this case (it was my understanding that using a Krylov-solver, only the matrixproduct Ax has to be provided and no matrix A needs to be assembled/saved in memory)?

Thanks

commented Sep 9, 2015 by maartent FEniCS User (3,910 points)

Btw I am stuck on a CentOS 6 system, hence version 1.4.0 (thanks to the conda packages of FEniCS). I noticed a new install option that uses HashDist: if there is any way to install 1.6.0 with this, I would be happy to hear so (gnu c library version: 2.12; gcc version:4.4.7)

commented Sep 9, 2015 by maartent FEniCS User (3,910 points)

You can try

mat = as_backend_type(A_B).mat()
mat.setUp()
mat.assemblyBegin()
mat.assemblyEnd()

Also try using petsc4py solver:

from petsc4py import PETSc
ksp = PETSc.KSP().create()
ksp.setOperators(mat, as_backend_type(P).mat())

PETSc.Options().setValue("ksp_type", "fgmres")
PETSc.Options().setValue("ksp_monitor", "")
PETSc.Options().setValue("ksp_converged_reason", "")
PETSc.Options().setValue("pc_type", "gamg")
ksp.setFromOptions()

x = as_backend_type(solution.vector()).vec()
b = as_backend_type(C).vec()

ksp.solve(b, x)

commented Sep 10, 2015 by Magne Nordaas FEniCS Expert (13,820 points)

Both work flawlessly, thanks (again)!

commented Sep 11, 2015 by maartent FEniCS User (3,910 points)

PETScKrylovSolver won't solve in version 1.4.0

1 Answer