[petsc-users] Implementing a homotopy solver

[zakaryah]

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).

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

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180716/f2583382/attachment.html>