# Diagnose "Newton Solver Did Not Converge" in FEniCS: Read the Iteration History First
Use this when a nonlinear solve dies with `Newton solver did not converge
because maximum number of iterations reached`. The iteration history —
absolute/relative residual per iteration — distinguishes a bad initial guess,
a wrong Jacobian, a scaling problem, and a genuinely hard nonlinearity, and
each needs a different fix. Covers legacy FEniCS (DOLFIN 2019.1.0,
`NonlinearVariationalSolver`) and FEniCSx (DOLFINx >= 0.6,
`dolfinx.nls.petsc.NewtonSolver`); the APIs differ and are labeled below.
## 1. Capture the error and turn on reporting
Legacy DOLFIN error:
```
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
```
DOLFINx error:
```
RuntimeError: Newton solver did not converge because maximum number of iterations reached
```
Before changing anything, make the solver print its residual history.
Legacy:
```python
prm = solver.parameters
prm["newton_solver"]["report"] = True
set_log_level(PROGRESS)
```
DOLFINx:
```python
from dolfinx.nls.petsc import NewtonSolver
solver = NewtonSolver(mesh.comm, problem)
solver.report = True
solver.convergence_criterion = "residual" # or "incremental"; see step 3
```
You will get one line per iteration. Legacy prints e.g.
`Newton iteration 0: r (abs) = 3.212e-01 (tol = 1.000e-10) r (rel) = 1.000e+00
(tol = 1.000e-09)`; DOLFINx prints e.g.
`Newton iteration 2: r (abs) = 5.00946 (tol = 1e-10) r (rel) = 0.285564 (tol = 1e-06)`.
## 2. Classify from the residual history
- **Residual decreases steadily, then stalls above tolerance (e.g. stuck at
1e-4 for many iterations): inconsistent Jacobian.** The search direction is
wrong, so Newton loses quadratic convergence. This is the most common cause
on the Discourse, and it almost always means a hand-derived Jacobian is
missing terms or is stale. Go to step 3.
- **Residual blows up or is NaN from iteration 0–1: bad initial guess or a
formulation bug.** The starting point is outside the basin of attraction,
or `F` itself is wrong (sign error, missing boundary term, unphysical
parameter). Go to step 4.
- **Residual oscillates or overshoots (down, up, down): needs damping.**
Go to step 5.
- **Residual reaches a small value but the solver still reports failure:
wrong convergence criterion or tolerances.** In mixed problems (e.g.
velocity–pressure), one field's scale dominates the residual norm, so the
"residual" criterion can never be satisfied even though the solution is
fine. Switch DOLFINx to `solver.convergence_criterion = "incremental"`
(norm of the update instead of the residual) and check whether the update
norm is genuinely small.
## 3. If the Jacobian is suspect: test it numerically
A Jacobian that is not the exact derivative of the residual destroys Newton
convergence. Verify with a finite-difference directional check (legacy
DOLFIN; `assemble`, `axpy`, `norm`, and matrix-vector `*` are all real
DOLFIN calls):
```python
import numpy as np
eps = 1e-7
r0 = assemble(F) # residual vector at current u
w = r0.copy(); w[:] = np.random.rand(w.local_size()) # random direction
J = assemble(derivative(F, u, TrialFunction(V))) # candidate Jacobian
u.vector().axpy(eps, w) # u -> u + eps*w
r1 = assemble(F)
u.vector().axpy(-eps, w) # restore u
fd = r1.copy(); fd.axpy(-1.0, r0) # r(u+eps*w) - r(u)
Jw = J * w
diff = fd.copy(); diff.axpy(-eps, Jw)
print("relative Jacobian error:", diff.norm("l2") / fd.norm("l2"))
```
A relative error near machine precision means the Jacobian is exact and the
problem lies elsewhere. A large error means the hand-written `J` is wrong —
replace it with the automatic derivative:
```python
# Legacy DOLFIN: let UFL differentiate instead of hand-coding J
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u_, bcs, J)
```
```python
# DOLFINx: NonlinearProblem builds J = derivative(F, u) internally
problem = fem.petsc.NonlinearProblem(F, uh, bcs=bcs)
```
Also check for non-smooth terms: `ufl.conditional`, `abs`, `max`/`min` have
zero or undefined UFL derivatives at kinks, and Newton stalls near them.
If your residual contains a conditional (contact, plasticity, phase-field),
expect stalling and consider smoothing or a linesearch (step 5).
## 4. If it blows up immediately: fix the starting point, not the solver
- **Warm-start from the linear problem.** Solve the linearized (or
small-load) version first and use it as the initial guess for the full
nonlinear problem. For Navier–Stokes, the standard warm start is the
Stokes solution.
- **Continuation / load stepping.** Ramp the hard parameter (inlet velocity,
pressure, Rayleigh number) in small steps, solving at each step and using
the previous solution as the guess. The Discourse Navier–Stokes threads
that "converge at 0.01 m/s but not 0.1 m/s" are almost always fixed by
continuation, not by solver tuning.
- **Check the formulation on a coarse mesh first.** If iteration 0 already
gives NaN, assemble `F` alone and inspect: `assemble(F).norm("l2")` at the
initial guess should be finite. NaN here means an assembly bug (division by
zero, `ln` of a non-positive quantity, uninitialized `Function`), not a
solver problem.
## 5. If it oscillates: damp the update
Reduce the step with the relaxation parameter (1.0 = full Newton step).
Legacy:
```python
prm["newton_solver"]["relaxation_parameter"] = 0.5
prm["newton_solver"]["maximum_iterations"] = 50
```
DOLFINx:
```python
solver.relaxation_parameter = 0.3 # < 1 damps; start at 0.5 and lower
solver.max_it = 500
```
Tune the tolerances deliberately instead of accepting defaults. Legacy:
`prm["newton_solver"]["absolute_tolerance"] = 1e-8`,
`prm["newton_solver"]["relative_tolerance"] = 1e-7`. DOLFINx:
`solver.atol`, `solver.rtol` (defaults are 1e-10 absolute / 1e-9 relative
scale — tighten or loosen based on the plateau value you saw in step 2).
## 6. Checklist
1. Enable reporting and capture the full residual history.
2. Classify: stalls (Jacobian) / blows up (initial guess or formulation) /
oscillates (damping) / tiny residual but "failure" (criterion/tolerance).
3. If Jacobian suspect: run the finite-difference check; use
`derivative(F, u, du)` instead of a hand-written `J`.
4. If initial guess: warm-start from the linear solution and/or add
continuation on the hard parameter.
5. If oscillation: lower `relaxation_parameter`, raise `maximum_iterations`
/ `max_it`, pick `residual` vs `incremental` to match the problem scaling.
6. Re-run with reporting on: a fixed problem shows the residual dropping
superlinearly in the last iterations — that signature is how you confirm
the fix rather than getting lucky.