[petsc-users] A number of questions about DMDA with SNES and Quasi-Newton methods

[Barry Smith]

> On Aug 26, 2017, at 5:56 PM, zakaryah . <zakaryah at gmail.com> wrote:
> 
> I'm using PETSc's SNES methods to solve PDEs which result from Euler-Lagrange equations for the strain energy of a 3D displacement field. There is an additional term in the Lagrangian which describes external forces which arise from various data sets, and that term contains nonlinearities (field terms higher than linear). The grid has about 1.6e6 elements, and the displacement field has 3 components at each grid element.
> 
> I'm trying to solve a sequence of successively more complicated equations, and the latest equation is failing to converge on some data sets. In particular, the methods were successful for the infinitesimal bulk strain (compression) energy, as well as the full infinitesimal strain energy (bulk + shear), but I'm now trying to generalize to the finite strain, as certain data sets are known to result from displacement fields for which the infinitesimal strain is a poor approximation.
> 
> I'm using a DMDA, closely following example 48, and my preferred solver is L-BFGS.

   So you are using ?

-snes_type qn -snes_qn_type lbfgs
 
> 
> I have read the FAQs "Why is Newton's method not converging?​" and "Why is my iterative linear solver not converging?​"​ which have raised a number of questions:

  Quasi Newton methods either don't use Jacobians or use only the initial Jacobian (the idea behind quasi-Newton methods is to approximate Jacobian information from previous iterations without having the user compute a Jacobian at each iteration). With PETSc's qn it only uses the Jacobian if you use the option

   -snes_qn_scale_type Jacobian

otherwise the Jacobian is never computed or used


> 
> Is there documentation for the DMDA/SNES methods somewhere?  I don't understand these very well.  For example, I am not allocating any matrix for the global Jacobian, and I believe this prevents me from changing the line search.  If I'm mistaken I would love to see an example of changing the line search type while using DMDA/SNES.

   Whether you provide a Jacobian or not is orthogonal to the line search. 

    You should be able to change the line search with 

-snes_linesearch_type bt  or nleqerr or basic or  l2  or cp

not all of them may work with qn


> 
> I don't know how to interpret the linesearch monitor.  Even for problems which are converging properly, the linesearch monitor reports "lssucceed=0" on every iteration.  Is this a problem?

  It returns a 0 if the line search does not believe it has achieved "sufficient decrease" in the function norm (or possibly some other measure of decrease) you should run -snes_linesearch_monitor also with the option -snes_monitor to see what is happening to the function norm

   For qn you can add the option

 -snes_qn_monitor

to get more detailed monitoring


> 
> I'm also having trouble understanding the methods for troubleshooting.  I suspect that I've made an error in the analytical Jacobian, which has a rather large number of non-zero elements, but I have no idea how to use -snes_type test -snes_test_display.  The FAQs mention that some troubleshooting tools are more useful for small test problems.  How small is small?

   Tiny, at most a few dozen rows and columns in the Jacobin.

   You should run without the -snes_test_display information, what does it say? Does it indicate the Jacobian or report there is likely a problem? 

    With DMDA you can also use -snes_fd_color to have PETSc compute the Jacobian for you instead of using your analytical form. If it works with this, but not your Jacobian then your Jacobian is wrong.

>  When I try to run the program with -snes_type test -snes_test_display, I get errors like:
> 
> [0]PETSC ERROR: Argument out of range [0]PETSC ERROR: Local index 1076396032 too large 4979879 (max) at 0
> 
> The second size is 1 less than the number of field elements, while the first number seems too large for any aspect of the problem - the Jacobian has at most 59 non-zero columns per row.
> 
> Because I suspect a possible error in the Jacobian, I ran with -snes_mf_operator -pc_type ksp -ksp_ksp_rtol 1e-12 and observed very similar failure to converge (diverging residual) as with the explicit Jacobian.

   What do you get with -ksp_monitor -ksp_ksp_monitor   it sounds like the true Jacobian is either very ill-conditioned or your function evaluation is wrong.

> Do I need to set an SNES method which is somehow compatible with the "matrix-free" approach? If I instead use -snes_mf, the problem seems to converge, but horrendously slowly (true residual relative decrease by about 1e-5 per iteration).  I suppose this supports my suspicion that the Jacobian is incorrect but doesn't really suggest a solution.
> 
> Is it possible that the analytical Jacobian is correct, but somehow pathological, which causes the SNES to diverge?

   Yes

>  Neither the Jacobian nor the function have singularities.
> 
> Thanks for any help you can provide!

    Try really hard to set up a small problem (like just use a very coarse grid) to experiment with as you try to get convergence. Using a big problem for debugging convergence is a recipe for pain.

   Also since you have a Jacobian I would start with -snes_fd_color -snes_type ls -snes_monitor -snes_linesearch_monitor -ksp_monitor -pc_type lu (not on a huge problem), what happens? Send the output 

   Barry