It appears to me that passing in boundary conditions in the solve command when using pre-assembled matrices is not working correctly. Consider the following problem: ************************************************** mesh = UnitSquareMesh(10,10) V = FunctionSpace(mesh,'CG',1) x,y = SpatialCoordinate(mesh) f_ex = 10.*exp(-(x-0.5)**2/0.02 - (y-0.5)**2/0.02) f = interpolate(f_ex,V) g_ex = sin(5*x) g = interpolate(g_ex,V) u, v = TrialFunction(V), TestFunction(V) bc = DirichletBC(V,0.0,(1,2)) A = assemble(inner(grad(u),grad(v))*dx,bcs=bc) b = assemble(f*v*dx + g*v*ds((3,4)),bcs=bc) U_1 = Function(V) solve(A,U_1,b) A = assemble(inner(grad(u),grad(v))*dx) b = assemble(f*v*dx + g*v*ds((3,4))) U_2 = Function(V) solve(A,U_2,b,bcs=bc) A = inner(grad(u),grad(v))*dx b = f*v*dx + g*v*ds((3,4)) U_3 = Function(V) solve(A==b,U_3,bcs=bc) **************************************************** The solve commands for U_1 and U_3 both converge, and the solutions agree to within machine precision. However, Firedrake fails to solve for U_2, with the error message: ConvergenceError: ('LinearSolver failed to converge after %d iterations with reason: %s', 1320, 'DIVERGED_DTOL' I suspect that boundary conditions are not being applied in the second solve command to the pre-assembled system. But according to the documentation on firedrakeproject.org, this should work as coded above. Am I doing something incorrectly? I'm also wondering how this would affect LinearVariationalSolver if there is indeed a problem with the boundary conditions. Thanks, Peter Sentz