[petsc-users] Fwd: Implementing a homotopy solver

[Matthew Knepley]

On Fri, Sep 28, 2018 at 11:13 PM zakaryah <zakaryah at gmail.com> wrote:

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

We need to clear up the usage first. If you want EXACTLY the same solver
for both solvers, then reuse
the KSP, otherwise do not do it. Does it work then?

  Thanks,

     Matt


> 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?
>>
>>

-- 
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.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180929/0f099e5b/attachment-0001.html>