[petsc-users] Implementing a homotopy solver

[Matthew Knepley]

On Mon, Jul 16, 2018 at 11:33 PM zakaryah <zakaryah at gmail.com> wrote:

> I'm still working on this, and I have a separate question about the
> implementation, which is related to calculating and removing the nullspace
> from the Jacobian.
>
> I have a system of n nonlinear equations in n+1 variables, due to the
> auxiliary variable (homotopy variable).  Given a guess, z, I'd like to
> calculate a corrector (Newton iteration) dz which solves
>
> J(z) dz = -F(z)
>
> where J is the n x (n+1) Jacobian of F.  Since J has a rank-one nullspace,
> I'd like to calculate that nullspace in each iteration of the SNES, then
> find the minimal norm solution dz.  The tricks for finding the nullspace
> without messing up the nice structure of the n x n submatrix of J are in
> Watson's paper: add a "dummy" row to J to form a square matrix A:
>
> A^T = [J^T e_k]
>
> and a trivial RHS:
>
> b^T = [0 ... 0 b_k]
>
> where the value of b_k is essentially arbitrary but can be chosen to
> orient the solution in the direction of increasing arc length.  A vector in
> the rank-one nullspace of J, u, solves
>
> A u = b              (1)
>
> and a specific solution for the Newton iterate solves
>
> A v = c              (2)
>
> where c^T = [-F^T c_{n+1}]
>
> and c_{n+1} is also an arbitrary constant.  The desired Newton update dz
> is just the Graham-Schmidt orthogonalization of v wrt u.
>
> From the previous discussion I think A can be implemented as a shell
> matrix as a straightforward way to efficiently solve the linear systems.  I
> have three questions about the implementation of this within the SNES:
>
> 1) What is the best way to "hook" a calculation of the nullspace of J, via
> eq. (1), within the SNES?  c is formed in FormFunction(), A is formed in
> FormJacobian(), but I need an intermediate step to form b and solve for u,
> then, I assume, set it as the nullspace of A (see question 3).
>

Just calculate the nullspace inside FormJacobian().


> 2) Is there a way to reuse the KSP for (1) to efficiently solve (2), given
> the method which answers question 1) ?
>

Not really, unless you use some type of factorization for the
preconditioner.


> 3) Can I set u to be in the nullspace of A?  It's really in the nullspace
> of the submatrix J, while the full matrix A has no nullspace, but the point
> is that I don't care about the last row of A - that's the whole trick of
> the homotopy.  I guess my question is whether anything in the solver will
> complain that Au != 0, because for all other purposes I want to project u
> out of the solution to (2) and it seems like treating it as a nullspace
> will achieve that.
>

Yes, I think it should be fine.

   Matt

-- 
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/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180717/5934b93e/attachment.html>