[petsc-users] Fwd: Implementing a homotopy solver

[zakaryah]

OK, thanks.

I'm using a composite DM, DM_packer.  To make a separate KSP, I did the
following in the FormJacobian() routine, after assembling the Jacobian
matrix A and the RHS for the tangent vector, b:

   - KSPCreate(PETSC_COMM_WORLD,&ksp)
   - KSPSetDM(ksp,DM_packer)
   - KSPSetDMActive(ksp,PETSC_FALSE) - because I want to set the operators
   - KSPSetOperators(ksp,A,P)
   - VecSet(n,0) - set initial guess to zero
   - KSPSolve(ksp,b,n) - this solve works correctly
   - VecNormalize(n,NULL)
   -
   MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp)
   - MatSetNullSpace(A,nullsp)

Then, with -snes_type newtonls -pc_type none -ksp_monitor
-ksp_monitor_true_residual -ksp_view, the output for the KSPSolve described
above, i.e. for the tangent vector, looks correct:

  0 KSP preconditioned resid norm 1.000000000000e+03 true resid norm
1.000000000000e+03 ||r(i)||/||b|| 1.000000000000e+00
...

185 KSP preconditioned resid norm 9.900713131874e-03 true resid norm
9.900713131904e-03 ||r(i)||/||b|| 9.900713131904e-06

Linear solve converged due to CONVERGED_RTOL iterations 185

KSP Object: 1 MPI processes

  type: gmres

    restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization
with no iterative refinement

    happy breakdown tolerance 1e-30

  maximum iterations=10000, initial guess is zero

  tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.

  left preconditioning

  using PRECONDITIONED norm type for convergence test

PC Object: 1 MPI processes

  type: none

  linear system matrix = precond matrix:

  Mat Object: 1 MPI processes

    type: seqaij

    rows=78247, cols=78247

    total: nonzeros=6063481, allocated nonzeros=6063481

    total number of mallocs used during MatSetValues calls =0

      using I-node routines: found 26083 nodes, limit used is 5


But the next KSP, i.e. the one in the SNES, doesn't converge in the true
residual:


  0 KSP preconditioned resid norm 2.045599092896e-04 true resid norm
2.803828296212e-04 ||r(i)||/||b|| 1.000000000000e+00

191 KSP preconditioned resid norm 2.009941278534e-09 true resid norm
1.636010142734e-04 ||r(i)||/||b|| 5.834915586465e-01


KSP Object: 1 MPI processes

  type: gmres

    restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization
with no iterative refinement

    happy breakdown tolerance 1e-30

  maximum iterations=10000, initial guess is zero

  tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.

  left preconditioning

  using PRECONDITIONED norm type for convergence test

PC Object: 1 MPI processes

  type: none

  linear system matrix = precond matrix:

  Mat Object: 1 MPI processes

    type: seqaij

    rows=78247, cols=78247

    total: nonzeros=6063481, allocated nonzeros=6063481

    total number of mallocs used during MatSetValues calls =0

      has attached null space

      using I-node routines: found 26083 nodes, limit used is 5


My guess about what's going on is that the tangent vector n isn't really in
the nullspace of A.  Rather, it's in the nullspace of the m x (m+1)
submatrix of A.  So, An=c e_{m+1}, where c is an arbitrary constant and
e_{m+1} is the m+1 th basis vector.  The nonlinear function also has an
added row, F_{m+1}, which is set to zero in the FormFunction() routine.  I
don't care about the value of F_{m+1}, but I suppose that if it's included
in the true residual, and NOT in the "preconditioned" residual, even with
pc_type none, then it will be tricky to diagnose the performance.  Should I
be using my own routine to evaluate the residual, so that F_{m+1} is not
included?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180930/02eb0615/attachment.html>