Hi, I'm working on a problem that requires solving a linear problem with non-constant bilinear and linear forms. I would like to wrap it into a LinearVariationalSolver, but also be able to access the assembled matrix and right hand side vector. I have figured out how to retrieve the assembled matrix, but as for grabbing the assembled right hand side, I've only been able to get it with a workaround. Here is a simple example of what I've been doing: ************************************************************************** from firedrake import * from firedrake.petsc import PETSc mesh = UnitSquareMesh(5,5) x,y = SpatialCoordinate(mesh) V = FunctionSpace(mesh,'CG',1) f = 1 + x*y u = TrialFunction(V) v = TestFunction(V) mu = Constant(0.0) k = Constant(0.0) a = inner(grad(u),grad(v))*dx + mu*u*v*dx L = k*f*v*dx U = Function(V) problem = LinearVariationalProblem(a,L, U,constant_jacobian=False) solver = LinearVariationalSolver(problem, solver_parameters={'ksp_type':'cg'}) # assign values to mu and k and solve mu.assign(1.0) k.assign(3.0) solver.solve() # grab Jacobian and check that it agrees with assembly of bilinear form A = solver.snes.getJacobian()[0].copy() A_ = assemble(a, mat_type='aij') A_.force_evaluation() A_ = A_.petscmat print((A - A_).norm()) # this prints 0.0, so they agree # grab rhs and check that it agrees with assembly of bilinear form # since F(u;v) = a(u,v) - L(v), and intial guess of function is zero, # then -F(0;v) = -a(0,v) + L(v) = L(v) b = -solver.snes.getFunction()[0].copy() b_ = assemble(L).dat.data b_ = PETSc.Vec().createWithArray(b_) print((b-b_).norm()) # this prints 6.5e-17, so they agree # change constants to new values and solve mu.assign(5.0) k.assign(-1.0) solver.solve() # grab Jacobian A = solver.snes.getJacobian()[0].copy() A_ = assemble(a, mat_type='aij') A_.force_evaluation() A_ = A_.petscmat print((A - A_).norm()) # prints 0.0, so they agree # grab rhs b = -solver.snes.getFunction()[0].copy() b_ = assemble(L).dat.data b_ = PETSc.Vec().createWithArray(b_) print((b-b_).norm()) # this prints 3.37 - they no longer agree # the problem is that the initial guess of the function is no longer zero, instead its the solution for the previous parameters # so, reset function U to have zero values U.assign(0) solver.solve() # check rhs again b = -solver.snes.getFunction()[0].copy() b_ = assemble(L).dat.data b_ = PETSc.Vec().createWithArray(b_) print((b-b_).norm()) # this prints 1.6e-17 - so they agree once more ************************************************************************************* So to grab the assembled vector, I need to make sure the solution is initialized to zero, and then evaluate the negative of the residual. Is there a cleaner way to grab the assembled vector out of the solver? Thanks, Peter Sentz