Why is Newton's method not converging?

Question

I am using PETSc's nonlinear solver package SNES to solve a system of nonlinear equations obtained by discretizing a partial differential equation. How can I determine why the solver is not converging and what can I do to successfully solve my equations?

Answer 1

Newton's method may not converge for many reasons, here are some of the most common.

• The Jacobian is wrong (or correct in sequential but not in parallel).
• The linear system is not solved or is not solved accurately enough.
• The Jacobian system has a singularity that the linear solver is not handling.
• There is a bug in the function evaluation routine.
• The function is not continuous or does not have continuous first derivatives (e.g. phase change or TVD limiters).
• The equations may not have a solution (e.g. limit cycle instead of a steady state) or there may be a "hill" between the initial guess and the steady state (e.g. reactants must ignite and burn before reaching a steady state, but the steady-state residual will be larger during combustion).

Here are some of the ways to help debug lack of convergence of Newton.

• Run with the options -snes_monitor -ksp_monitor_true_residual -snes_converged_reason -ksp_converged_reason. If the linear solve does not converge, check if the Jacobian is correct, then see this question. If the preconditioned residual converges, but the true residual does not, the preconditioner may be singular. If the linear solve converges well, but the line search fails, the Jacobian may be incorrect.
• Run with -pc_type lu or -pc_type svd to see if the problem is a poor linear solver
• Run with -mat_view or -mat_view_draw to see if the Jacobian looks reasonable
• Run with -snes_type test -snes_test_display to see if the Jacobian you are using is wrong. Compare the output when you add -mat_fd_type ds to see if the result is sensitive to the choice of differencing parameter.
• Run with -snes_mf_operator -pc_type lu to see if the Jacobian you are using is wrong. If the problem is too large for a direct solve, try -snes_mf_operator -pc_type ksp -ksp_ksp_rtol 1e-12. Compare the output when you add -mat_mffd_type ds to see if the result is sensitive to choice of differencing parameter.
• Run on one processor to see if the problem is only in parallel.
• Run with -snes_ls_monitor to see if the line search is failing (this is usually a sign of a bad Jacobian).
• Run with -info to get more detailed information on the solution process.

Here are some ways to help the Newton process if everything above checks out

• Run with grid sequencing (-snes_grid_sequence is all you need if working with a DM) to generate better initial guess on your finer mesh
• Run with quad precision (./configure --with-precision=__float128 --download-f2cblaslapack with PETSc 3.2 and later, needs version 4.6 or later GNU compilers)
• Change the units (nondimensionalization), boundary condition scaling, or formulation so that the Jacobian is better conditioned.
• Mollify features in the function that do not have continuous first derivatives (often occurs when there are if statements in the residual evaluation, e.g. phase change or TVD limiters). Use a variational inequality solver (SNESVINEWTONRSLS) if the discontinuities are of fundamental importance.
• Try a trust region method (-ts_type tr, may have to adjust parameters).
• Run with some continuation parameter from a point where you know the solution, see TSPSEUDO for solving steady-state problems. There are homotopy solver packages like PHCpack that can get you all possible solutions (and tell you that it has found them all) but those are not scalable and cannot solve anything but small problems.