[petsc-users] Fwd: Implementing a homotopy solver

[zakaryah]

I'm working on a homotopy solver which follows a zero curve through
solution space.  A central aspect of this method is to calculate the vector
tangent to the curve, predict the next point, then correct using iteration
of, e.g. Newton's method, orthogonal to the tangent vector.

Previously, we discussed the possibilities of implementing this within
PETSc's SNES.  Within my FormJacobian() function, I construct the linear
system which defines the tangent vector, solve it, then add the vector to
the nullspace of the Jacobian.  I think that in principle this can work,
but I suspect I'm doing something wrong.

Here's a summary of the code within FormJacobian():

   - Set values and assemble Jacobian matrix A - this is working fine
   - Set values and assemble RHS vector b for linear system defining
   tangent vector n - this is working fine
   - SNESGetKSP(snes,&my_ksp) - I thought it made sense to use the KSP
   associated with the SNES, hoping that PCs which use a factorization could
   be reused when the SNES calls KSPSolve() to calculate the update
   - KSPSetFromOptions(my_ksp) - not sure this is necessary but one of my
   problems is setting options for this KSP from the command line and even
   with this call it doesn't seem to be working properly
   - MatSetNullSpace(A,NULL) - remove any existing null space from Jacobian
   - KSPSetOperators(my_ksp,A,P) - P is the other matrix in FormJacobian()
   - VecSet(n,0) - set initial guess to zero
   - KSPSolve(my_ksp,b,n) - these solves appear to work, i.e. use the
   options passed from the command line with -ksp_XXX or -pc_XXX
   - VecNormalize(n,NULL)
   -
   MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp)
   - MatSetNullSpace(A,nullsp)
   - MatNullSpaceDestroy(&nullsp)
   - return

The immediate problem is that the subsequent KSPSolve(), i.e. the one
called internally by SNESSolve(), behaves strangely.  For example, if I use
-pc_type none -ksp_monitor -ksp_monitor_true_residual, then the KSPSolve()
that I call within FormJacobian() looks correct - "preconditioned" norm and
true norm are identical, and both converge as I expect (i.e. slowly but
geometrically).  However, the subsequent KSPSolve(), internal to the
SNESSolve(), has large differences between the preconditioned norm and the
true norm.  In addition, the KSP does not converge in the true residual,
but I'll have a hard time debugging that without knowing how to properly
set the options.

I hope someone can help me see what I'm doing wrong.

On Sun, Jul 22, 2018 at 9:09 PM zakaryah <zakaryah at gmail.com> wrote:

> ​Thanks Matt and Barry,
>
> Matt - if I do the calculation in FormJacobian(), which makes by far the
> most sense and is as per your suggestion, do I need to set the operators of
> the SNES's KSP back to whatever they were before I set them?  The linear
> system I want to solve within FormJacobian() involves the Jacobian matrix
> itself, and I want to remove the "nullspace" from that same matrix within
> FormFunction().
>
> Barry - I'm trying to implement a homotopy solver.  In short, I have a
> system of n nonlinear equations in n variables, F(x), which is hard to
> solve because the Jacobian tends to become singular using standard
> methods​.  I want to add an auxiliary variable, lambda, to create a
> homotopy: H(lambda,x) = lambda*F(x) + (1-lambda)G(x), where G is "easy to
> solve", and the idea is that the Jacobian of the n+1 variable system will
> not become singular along the curve H(lambda,x) = 0.
>
> The method involves adding an equation to H so that the Jacobian H' is
> square.  The "submatrix" refers to the n x (n+1) matrix which represents
> the Jacobian without the added equation, whereas my FormJacobian() forms
> the entire (n+1) x (n+1) matrix H'.  I only refer to the submatrix because
> it has a nullspace, and I want to find it by solving a linear system
> designed for this purpose, H' u = b, where b is not the zero vector.  H'
> has no nullspace, but I want to remove the projection of u from my SNES
> solution vector, as u IS in the nullspace of the submatrix.
>
> The RHS vector b can be calculated outside the SNES solve.  I guess my
> FormJacobian() should look like this:
>
> FormJacobian(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx) {
>   KSP ksp;
>   MatNullSpace unull;
>   user_struct *user = (user_struct*)ctx;
>
>   calculate Amat
>   assemble Amat
>
>   if (Pmat != Amat) {
>     assemble Amat
>   }
>
>   SNESGetKSP(snes,&ksp);
>   KSPSetOperators(ksp,Amat,Pmat);
>   KSPSolve(ksp,user->b,user->u);
>
>   MatNullSpaceCreate(PetscObjectComm((PetscObject)Amat), PETSC_FALSE, 1,
> &(user->u),&unull);
>   MatSetNullSpace(Amat,unull);
>   MatNullSpaceDestroy(&unull);
> }
>
> Does this look right?
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180928/c10077c3/attachment.html>