Solving pre-assembled linear system: adjoint not working

[nbouziani]

Despite having a dedicated block and annotation mechanism, it seems that solving a pre-assembled linear system (i.e. via solve(A, x, b)) is not adjoint compatible. The corresponding block seems wrong, e.g it assumes the rhs is not a ufl.Coefficient (cf. here). Also, the pyadjoint core block class for solve (i.e. GenericSolveBlock), which is currently used to compute the adjoint/tlm of solving pre-assembled systems, is not equipped for that case.

I think this has been going off the radar because there are no tests hitting that implementation.

Example:

from firedrake import *
from firedrake_adjoint import *

mesh = IntervalMesh(10, 0, 1)
V = FunctionSpace(mesh, "Lagrange", 1)

f = Function(V)
f.vector()[:] = 1

u = TrialFunction(V)
v = TestFunction(V)
bc = DirichletBC(V, Constant(1), "on_boundary")

a = inner(grad(u), grad(v))*dx
A = assemble(a, bcs=bc)

L = f * v * dx
b = assemble(L)

u_ = Function(V)
# This doesn't work (relies on `SolveLinearSystemBlock`)
solve(A, u_, b)
# This works (relies on `SolveVarFormBlock`)
# solve(a == L, u_, bcs=bc)

J = assemble(u_**2*dx)
Jhat = ReducedFunctional(J, Control(f))
Jhat.derivative()