Dear All, I have trilinear form a(w,u,v) and usually solve a(w,u,v) = l(v), for u with fixed w. Here I have: a = w dot(grad(u), grad(v)) * dx, and l(v) = f * v * dx. Now, given u_tilde I want to solve for w, i.e. find w_tilde such, that a(w_tilde, u_tilde, v) = l(v), for any v in suitable subspace. I'm trying the following: mesh = UnitSquareMesh(size, size) V = FunctionSpace(mesh, "CG", 1) f = Function(V) f.interpolate(Expression("1")) u_tilde = Function(V) u_tilde.interpolate(Expression("1/(1+exp(-1*x[0])) * 1/(1+exp(-1*(1-x[0]))) *1/(1+exp(-1*x[1])) * 1/(1+exp(-1*(1-x[1])))" )) w_tilde = TrialFunction(V) v = TestFunction(V) a_tilde = (w_tilde * dot(grad(u_tilde), grad(v))) * dx frm_tilde = assemble( a_tilde ) L = f * v * dx rhs = assemble(L) w_sol = Function(V) solve(frm_tilde, rhs, w_sol, solver_parameters={"ksp_type": "preonly", "pc_type": "lu"}) File("w_sol.pvd") << w_sol It fails to converge: 140 r = self.ksp.getConvergedReason() 141 if r < 0: --> 142 raise RuntimeError("LinearSolver failed to converge after %d iterations with reason: %s", self.ksp.getIterationNumber(), solving_utils.KSPReasons[r]) KeyError: -11 I'm using direct method, not the iterative one, so what is happening here?