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?
-
8$\begingroup$ As with any iterative method, it is paramount that one come up with a good seed to start up the Newton-Raphson method. A poor starting point often results in chaos. $\endgroup$– J. M.Nov 30, 2011 at 7:57
-
6$\begingroup$ I think that "often results in chaos" is incorrect. Referring to Blum, Cucker, Shub, and Smale, Newton has attractive basins separated by boundaries which may result in chaotic iteration. Thus this kind of behavior is very unlikely compared to convergence. If the algorithm is only looking for real solutions, it will fail to converge often, but will not become chaotic. $\endgroup$– Matt KnepleyDec 1, 2011 at 13:30
-
6$\begingroup$ I meant "chaos" in the nontechnical sense, @Matt; for instance, divergence to infinity, or cyclic behavior. Maybe "disaster" would have been a better choice of words? $\endgroup$– J. M.Dec 2, 2011 at 6:45
1 Answer
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 aDM
) 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.