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

[{"constraint": "Any legacy dolfin release (e.g. 2019.1.0); PETScKrylovSolver parameters absolutetolerance, relativetolerance, maximumiterations", "technology": "Legacy FEniCS (DOLFIN) PETScKrylovSolver", "versionscheme": "unknown"}, {"constraint": ">=0.6.0; applies when an iterative KSP is selected via petscoptions or built manually (LinearProblem defaults to direct LU)", "technology": "FEniCSx (DOLFINx) PETSc KSP", "versionscheme": "semver"}]

Diagnose PETSc Krylov Solver Divergence in FEniCS: Read the PETSc Reason First

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

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