Diagnose Newton solver non-convergence in FEniCS: read the iteration history first

TL;DR

For "RuntimeError: Newton solver did not converge because maximum number of iterations reached": enable reporting and capture the full residual history. 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.

RuntimeError: Newton solver did not converge because maximum number of iterations reached

When to use

You are seeing this: Covers legacy FEniCS (DOLFIN 2019.1.0, NonlinearVariationalSolver) and FEniCSx (DOLFINx >= 0.6, dolfinx.nls.petsc.NewtonSolver); the APIs differ and are labeled below. Use this skill when you run into "Diagnose Newton solver non-convergence in FEniCS: read the iteration history first".

When not to use

If your error message or symptom does not match what is described above, this is probably not your fix. Search for your exact error text instead of forcing this one to fit.

Versions

Versions mentioned in the source: DOLFIN 2019.1.0. If you are on something much newer or older, the details may have shifted.

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:

prm = solver.parameters
prm["newton_solver"]["report"] = True
set_log_level(PROGRESS)

DOLFINx:

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):

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:

# Legacy DOLFIN: let UFL differentiate instead of hand-coding J
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u_, bcs, J)

# 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:

prm["newton_solver"]["relaxation_parameter"] = 0.5
prm["newton_solver"]["maximum_iterations"] = 50

DOLFINx:

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.

Variants

Other phrasings of the same problem that show up in reports:

*** Error: Unable to solve nonlinear system with NewtonSolver.

*** 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.

import numpy as np

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"))

Built from FEniCS Discourse threads on Newton divergence (Navier-Stokes at higher inlet velocity, hyperelasticity with mixed elements, manual Newton implementations) and the FEniCS tutorial's nonlinear-problem parameter documentation. Encodes the real solver-reporting APIs (legacy prm['newtonsolver'], DOLFINx solver.report/relaxationparameter/maxit/convergencecriterion) and the real failure message 'Newton solver did not converge because maximum number of iterations reached'.