Hello, I am a FEniCS user but new to petsc4py. I am trying to modify the KSP solver through the SNES object to implement the Sherman-Morrison formula (e.g. http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html ). I am solving a nonlinear system of the form [K(u)−kaaT]Δu=−F(u). Here the jacobian matrix K is modified by the term kaaT, where k is a scalar. Notably, K is sparse, while the term kaaT results in a full matrix. This problem can be solved efficiently using the Sherman-Morrison formula :
[K−kaaT]-1 = K-1 - (kK-1 aaTK-1)/(1+kaTK-1a) I have managed to successfully implement this at the linear solve level (by modifying the KSP solver) inside a custom Newton solver in python by following an incomplete tutorial at https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner : * while (norm(delU) > alpha): # while not converged * * self.update_F() # call to method to update r.h.s form * self.update_K() # call to update the jacobian form * K = assemble(self.K) # assemble the jacobian matrix * F = assemble(self.F) # assemble the r.h.s vector * a = assemble(self.a_form) # assemble the a_form (see Sherman Morrison formula) * * for bc in self.mem.bc: # apply boundary conditions * bc.apply(K, F) * bc.apply(K, a) * * B = PETSc.Mat().create() * * # Assemble the bilinear form that defines A and get the concrete * # PETSc matrix * A = as_backend_type(K).mat() # get the PETSc objects for K and a * u = as_backend_type(a).vec() * * # Build the matrix "context" # see firedrake docs * Bctx = MatrixFreeB(A, u, u, self.k) * * # Set up B * # B is the same size as A * B.setSizes(*A.getSizes()) * B.setType(B.Type.PYTHON) * B.setPythonContext(Bctx) * B.setUp() * * * ksp = PETSc.KSP().create() # create the KSP linear solver object * ksp.setOperators(B) * ksp.setUp() * pc = ksp.pc * pc.setType(pc.Type.PYTHON) * pc.setPythonContext(MatrixFreePC()) * ksp.setFromOptions() * * solution = delU # the incremental displacement at this iteration * * b = as_backend_type(-F).vec() * delu = solution.vector().vec() * * ksp.solve(b, delu) * * self.mem.u.vector().axpy(0.25, self.delU.vector()) # poor man's linesearch * counter += 1 Here is the corresponding petsc4py code adapted from the firedrake docs: 1. class MatrixFreePC(object): 2. 3. def setUp(self, pc): 4. B, P = pc.getOperators() 5. # extract the MatrixFreeB object from B 6. ctx = B.getPythonContext() 7. self.A = ctx.A 8. self.u = ctx.u 9. self.v = ctx.v 10. self.k = ctx.k 11. # Here we build the PC object that uses the concrete, 12. # assembled matrix A. We will use this to apply the action 13. # of A^{-1} 14. self.pc = PETSc.PC().create() 15. self.pc.setOptionsPrefix("mf_") 16. self.pc.setOperators(self.A) 17. self.pc.setFromOptions() 18. # Since u and v do not change, we can build the denominator 19. # and the action of A^{-1} on u only once, in the setup 20. # phase. 21. tmp = self.A.createVecLeft() 22. self.pc.apply(self.u, tmp) 23. self._Ainvu = tmp 24. self._denom = 1 + self.k*self.v.dot(self._Ainvu) 25. 26. def apply(self, pc, x, y): 27. # y <- A^{-1}x 28. self.pc.apply(x, y) 29. # alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u) 30. alpha = (self.k*self.v.dot(y)) / self._denom 31. # y <- y - alpha * A^{-1}u 32. y.axpy(-alpha, self._Ainvu) 33. 34. 35. class MatrixFreeB(object): 36. 37. def __init__(self, A, u, v, k): 38. self.A = A 39. self.u = u 40. self.v = v 41. self.k = k 42. 43. def mult(self, mat, x, y): 44. # y <- A x 45. self.A.mult(x, y) 46. 47. # alpha <- v^T x 48. alpha = self.v.dot(x) 49. 50. # y <- y + alpha*u 51. y.axpy(alpha, self.u) However, this approach is not efficient as it requires many iterations due to the Newton step being fixed, so I would like to implement it using SNES and use line search. Unfortunately, I have not been able to find any documentation/tutorial on how to do so. Provided I have the FEniCS forms for F, K, and a, I'd like to do something along the lines of: solver = PETScSNESSolver() # the FEniCS SNES wrapper snes = solver.snes() # the petsc4py SNES object ## ?? ksp = snes.getKSP() # set ksp option similar to above solver.solve() I would be very grateful if anyone could could help or point me to a reference or demo that does something similar (or maybe a completely different way of solving the problem!). Many thanks in advance! Alex
