https://gitlab.com/petsc/petsc/issues/557
> On Feb 5, 2020, at 7:35 AM, Matthew Knepley <[email protected]> wrote: > > Perhaps Barry is right that you want Picard, but suppose you really want > Newton. > > "This problem can be solved efficiently using the Sherman-Morrison formula" > Well, maybe. The main assumption here is that inverting K is cheap. I see two > things you can do in a straightforward way: > > 1) Use MatCreateLRC() > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateLRC.html > to create the Jacobian > and solve using an iterative method. If you pass just K was the > preconditioning matrix, you can use common PCs. > > 2) We only implemented MatMult() for LRC, but you could stick your SMW code > in for MatSolve_LRC if you really want to factor K. We would > of course help you do this. > > Thanks, > > Matt > > On Wed, Feb 5, 2020 at 1:36 AM Smith, Barry F. via petsc-users > <[email protected]> wrote: > > 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 <[email protected]> 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 > > > > -- > What most experimenters take for granted before they begin their experiments > is infinitely more interesting than any results to which their experiments > lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/
