Diagnose PETSc Krylov solver divergence in FEniCS: read the PETSc reason first

TL;DR

For "DIVERGEDPCFAILED": check the tolerances first, the cheapest fix. A real Discourse case set absolute\tolerance = 1e-17 with maximum\iterations = 10000 and got DIVERGED\_ITS with ||r|| = 1.365103e-11: the solver had converged to a perfectly good solution and only "failed" because the tolerance was tighter than floating point allows.

DIVERGED_PC_FAILED

When to use

You are seeing this: The parenthesized PETSc reason in the error, DIVERGED\ITS, DIVERGED\PC\FAILED, DIVERGED\NANORINF, tells you which of three unrelated problems you have. Use this skill when you run into "Diagnose PETSc Krylov solver divergence in FEniCS: read the PETSc reason 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: 2019.1.0. If you are on something much newer or older, the details may have shifted.

Use this when a linear solve fails with Unable to solve linear system using PETSc Krylov solver (legacy DOLFIN) or a PETSc KSP divergence (DOLFINx). The parenthesized PETSc reason in the error, DIVERGED_ITS, DIVERGED_PC_FAILED, DIVERGED_NANORINF, tells you which of three unrelated problems you have. Covers legacy DOLFIN (2019.1.0) and DOLFINx (>= 0.6). Note: DOLFINx's LinearProblem defaults to a direct LU solve (ksp_type=preonly, pc_type=lu) unless you override petsc_options, so in DOLFINx this skill applies when you chose an iterative solver yourself or built the KSP manually; in legacy DOLFIN iterative Krylov solvers are the common path.

1. Read the PETSc reason in the error

Legacy DOLFIN prints, for example:

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 1.365103e-11).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.

In DOLFINx, get the same information by asking the KSP for its converged reason after the solve, or set ksp_monitor/ksp_view in petsc_options:

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={
    "ksp_type": "gmres", "pc_type": "ilu",
    "ksp_monitor": None, "ksp_converged_reason": None})

The three reasons and what each means:

• DIVERGED_ITS, hit maximum_iterations with the residual still above tolerance. The solver ran fine but couldn't get there: wrong solver for the matrix, singular system, or absurd tolerances. Go to step 2.

• DIVERGED_PC_FAILED, failed in 0 iterations: the preconditioner setup failed (factorization broke down or ran out of memory). Go to step 3.

• DIVERGED_NANORINF, NaN or Inf in the matrix or right-hand side. This is an assembly bug, not a solver problem. Go to step 4.

2. DIVERGED_ITS: classify why the iterations couldn't converge

2a. Check the tolerances first, the cheapest fix. A real Discourse case set absolute_tolerance = 1e-17 with maximum_iterations = 10000 and got DIVERGED_ITS with ||r|| = 1.365103e-11: the solver had converged to a perfectly good solution and only "failed" because the tolerance was tighter than floating point allows. If your final residual norm is small (near 1e-10 or below), loosen the tolerances instead of changing the solver. Legacy:

solver = PETScKrylovSolver("cg", "ilu")
solver.parameters["absolute_tolerance"] = 1e-10
solver.parameters["relative_tolerance"] = 1e-6
solver.parameters["maximum_iterations"] = 1000

2b. Match the solver to the matrix. CG requires a symmetric positive-definite matrix. Advection terms, non-symmetric stabilization, or saddle-point structure make CG stall or diverge, switch to GMRES:

# Legacy DOLFIN
solver = PETScKrylovSolver("gmres", "ilu")
# or via options database
PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "ilu")

# DOLFINx
problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={
    "ksp_type": "gmres", "pc_type": "ilu"})

For symmetric indefinite saddle-point systems (Stokes, mixed Poisson), monolithic ILU preconditioning usually fails; use a block (fieldsplit) preconditioner or MINRES with a Schur-complement preconditioner instead of fighting the iteration count.

2c. Check for a singular matrix, pure Neumann or unconstrained rigid-body modes. If every boundary condition is Neumann (pressure Poisson in a projection method is the classic case), the matrix has a constant nullspace and Krylov solvers stall or return garbage. Attach the nullspace and orthogonalize the right-hand side. DOLFINx:

from petsc4py import PETSc
A = fem.petsc.assemble_matrix(fem.form(a), bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(fem.form(L))
fem.petsc.apply_lifting(b, [fem.form(a)], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, bcs)
nullspace = PETSc.NullSpace().create(constant=True, comm=mesh.comm)
assert nullspace.test(A)   # confirms A really has the constant nullspace
A.setNullSpace(nullspace)
nullspace.remove(b)        # b must have no nullspace component
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("cg")
ksp.getPC().setType("hypre")
ksp.setFromOptions()
uh = fem.Function(V)
ksp.solve(b, uh.vector)

Symptom that points here: the residual decreases for a while then plateaus far above tolerance, on a problem with no Dirichlet conditions (or with floating subdomains in elasticity, where the nullspace is the six rigid-body modes). Fixing the physics (adding a Dirichlet condition or a Lagrange multiplier pinning the mean) is equivalent and often cleaner.

2d. Check mesh-dependent blowup. If the same code converges on a coarse mesh and diverges on a fine one, suspect the preconditioner (ILU degrades with problem size) rather than the discretization, move to AMG (gamg / hypre) for elliptic problems before assuming the formulation is wrong.

3. DIVERGEDPCFAILED: the preconditioner broke, not the iteration

Zero iterations means factorization failed during setup. Common causes:

• Direct solver ran out of memory on a large 3D problem (MUMPS via pc_type=lu). The turtleFSI known-issues page documents exactly this failure; the remedy is to give MUMPS more workspace: PETScOptions.set("mat_mumps_icntl_14", 400) (percentage memory increase), or switch to an iterative solver with AMG.

• Factorization of a singular or indefinite matrix. ILU/LU on a matrix with a zero pivot (see 2c) fails at setup. Fix the singularity first.

• Preconditioner incompatible with the matrix type, e.g. ICC on a non-SPD matrix. Match pc_type to the solver choice from 2b.

4. DIVERGED_NANORINF: find the NaN in the assembly

The solver is innocent; something in a or L evaluated to NaN/Inf. Bisect: assemble L alone and check b.norm() is finite; then assemble each term of a as its own form and check norms. Typical sources: division by a quantity that is zero at some quadrature points (1/J with inverted elements, check mesh quality), ln of a non-positive argument, uninitialized Functions used as coefficients, or an Expression with a C-level domain error. Fix the term; no solver option helps here.

5. Checklist

1. Copy the PETSc reason from the error: DIVERGED_ITS / DIVERGED_PC_FAILED / DIVERGED_NANORINF.

2. DIVERGED_ITS: tolerances sane? (2a) solver matches matrix symmetry? (2b) singular system needing a nullspace? (2c) fine-mesh-only → preconditioner scaling? (2d).

3. DIVERGED_PC_FAILED: direct-solver memory (mat_mumps_icntl_14) or a singular/indefinite matrix breaking the factorization.

4. DIVERGED_NANORINF: bisect the assembly term by term; fix the NaN source.

5. Confirm the fix by re-running with ksp_monitor / ksp_converged_reason visible: the reason should read CONVERGED_RTOL (or CONVERGED_ATOL), not merely "no error raised".

Variants

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

*** Error: Unable to solve linear system using PETSc Krylov solver.

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 1.365103e-11).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.

Built from FEniCS Discourse threads on Krylov divergence (DIVERGED_ITS with over-tight absolute tolerance, pure-Neumann singular Poisson nullspace handling, MUMPS memory failures documented in turtleFSI known issues) and the real legacy DOLFIN error format 'Solution failed to converge in N iterations (PETSc reason ..., residual norm ||r|| = ...)'. Nullspace recipe verified against Discourse DOLFINx examples (PETSc.NullSpace().create(constant=True), setNullSpace, nullspace.remove(b)).