[petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?

[Smith, Barry F.]

   I am not sure of everything in your email but it sounds like you want to use a "Picard" iteration to solve [K(u)−kaaT]Δu=−F(u). That is solve

  A(u^{n}) (u^{n+1} - u^{n}) = F(u^{n}) - A(u^{n})u^{n}  where A(u) = K(u) - kaaT

 PETSc provides code to this with SNESSetPicard() (see the manual pages) I don't know if Petsc4py has bindings for this.

  Adding missing python bindings is not terribly difficult and you may be able to do it yourself if this is the approach you want.

   Barry

   

> On Feb 4, 2020, at 5:07 PM, Olek Niewiarowski <aan2 at princeton.edu> wrote:
> 
> 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:
> 
> 	• class MatrixFreePC(object):  
> 	•   
> 	•     def setUp(self, pc):  
> 	•         B, P = pc.getOperators()  
> 	•         # extract the MatrixFreeB object from B  
> 	•         ctx = B.getPythonContext()  
> 	•         self.A = ctx.A  
> 	•         self.u = ctx.u  
> 	•         self.v = ctx.v  
> 	•         self.k = ctx.k  
> 	•         # Here we build the PC object that uses the concrete,  
> 	•         # assembled matrix A.  We will use this to apply the action  
> 	•         # of A^{-1}  
> 	•         self.pc = PETSc.PC().create()  
> 	•         self.pc.setOptionsPrefix("mf_")  
> 	•         self.pc.setOperators(self.A)  
> 	•         self.pc.setFromOptions()  
> 	•         # Since u and v do not change, we can build the denominator  
> 	•         # and the action of A^{-1} on u only once, in the setup  
> 	•         # phase.  
> 	•         tmp = self.A.createVecLeft()  
> 	•         self.pc.apply(self.u, tmp)  
> 	•         self._Ainvu = tmp  
> 	•         self._denom = 1 + self.k*self.v.dot(self._Ainvu)  
> 	•   
> 	•     def apply(self, pc, x, y):  
> 	•         # y <- A^{-1}x  
> 	•         self.pc.apply(x, y)  
> 	•         # alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u)  
> 	•         alpha = (self.k*self.v.dot(y)) / self._denom  
> 	•         # y <- y - alpha * A^{-1}u  
> 	•         y.axpy(-alpha, self._Ainvu)  
> 	•   
> 	•   
> 	• class MatrixFreeB(object):  
> 	•   
> 	•     def __init__(self, A, u, v, k):  
> 	•         self.A = A  
> 	•         self.u = u  
> 	•         self.v = v  
> 	•         self.k = k  
> 	•   
> 	•     def mult(self, mat, x, y):  
> 	•         # y <- A x  
> 	•         self.A.mult(x, y)  
> 	•   
> 	•         # alpha <- v^T x  
> 	•         alpha = self.v.dot(x)  
> 	•   
> 	•         # y <- y + alpha*u  
> 	•         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