[petsc-users] Fwd: Implementing a homotopy solver

[zakaryah]

​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/20180722/3ddfd2fa/attachment.html>