Odd stuff with evaluation order and solvers

[Former user]

Hi!

I am trying to set up a petsc amg solver for 3D elasticity and ran into some problems with dolfin-adjoint giving an error. I have boiled it down to a minimal example with a Poisson equation:

from dolfin import *
from dolfin_adjoint import *

parameters["linear_algebra_backend"] = "PETSc"

mesh = UnitSquareMesh(8, 8)

V = FunctionSpace(mesh, 'Lagrange', 1)
Vd = FunctionSpace(mesh, 'DG', 0)
rho = interpolate(Constant(0.5), Vd, name="Density")

def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

u = TrialFunction(V)
v = TestFunction(V)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)
a = rho*inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
A, b = assemble_system(a, L, bc)

sol=Function(V)

# Setting solver does not work
pc = PETScPreconditioner("petsc_amg")
solver = PETScKrylovSolver("gmres", pc)
solver.set_operator(A)
solver.solve(sol.vector(), b)
#solve(A,sol.vector(),b,"gmres","petsc_amg") # Uncomment this line and the error is gone!

m=Control(rho)

#print compute_gradient(Functional(0.5*inner(sol,sol)*dx),m,forget=False) # Uncomment this line and the error is gone!
print compute_gradient(Functional(rho*dx),m,forget=False)
print compute_gradient(Functional(0.5*inner(sol,sol)*dx),m,forget=False)

The error I am getting is:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [e-mail]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to unable to solve linear system with PETSc Krylov solver.
*** Reason:  Non-matching dimensions for linear system (matrix has 0 rows and right-hand side vector has 81 rows).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2017.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

If the code is run with 'solve' and not with the 'solver class' it does not throw an error. Also if the last functional is evaluated first the code runs, even with the 'solver class'.

What am I doing wrong?

Best Regards, Søren

PS Thanks for some great software :-)

[Søren Madsen]

Hi again!

It seems like in the file 'petsc_krylov_solver.py' in the 'solve' method of the 'PETScKrylovSolver' class the matrices are not set correctly. Inserting the following print just before the solve prints out '(None,None)' for the second call to 'compute_gradient'.

                    print(solver.operators)
                    solver.solve(x.vector(), rhs)

My best guess is that this part

                       if adj_petsc_krylov_solvers[idx] is None:
                            need_to_set_operator = True
                            adj_petsc_krylov_solvers[idx] = PETScKrylovSolver(*solver_parameters)
                            adj_ksp = adj_petsc_krylov_solvers[idx].ksp()
                            fwd_ksp = petsc_krylov_solvers[idx].ksp()
                            adj_ksp.setOptionsPrefix(fwd_ksp.getOptionsPrefix())
                            adj_ksp.setType(fwd_ksp.getType())
                            adj_ksp.pc.setType(fwd_ksp.pc.getType())
                            adj_ksp.setFromOptions()
                        else:
                            need_to_set_operator = self._need_to_reset_operator

is not doing its job correctly as in the second call to 'compute_gradient' it goes into the 'else' part but the operator have never been set correctly for the second call to 'compute_gradient'. Will it work to just set the operator every time? I.e. change the above to

                       #if adj_petsc_krylov_solvers[idx] is None:
                       need_to_set_operator = True
                       adj_petsc_krylov_solvers[idx] = PETScKrylovSolver(*solver_parameters)
                       adj_ksp = adj_petsc_krylov_solvers[idx].ksp()
                       fwd_ksp = petsc_krylov_solvers[idx].ksp()
                       adj_ksp.setOptionsPrefix(fwd_ksp.getOptionsPrefix())
                       adj_ksp.setType(fwd_ksp.getType())
                       adj_ksp.pc.setType(fwd_ksp.pc.getType())
                       adj_ksp.setFromOptions()
                       #else:
                       #    need_to_set_operator = self._need_to_reset_operator

The above code shows no errors when doing this, but is it actually correct? :-)