Hi Peter, It appears our documentation is lagging development here. The ability to add boundary conditions after assembling a matrix was removed recently, basically because it causes a lot of complication under the hood and doesn’t really add any functionality. If you want to pre-assemble systems for some reason (and there is usually no reason to do this), then you should apply the boundary conditions at assembly time. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Sentz, Peter" <psentz@sandia.gov> Date: Monday, 1 July 2019 at 23:44 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Solve no longer applies BC's to pre-assembled matrices 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