# 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`:
```python
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:
```python
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:
```python
# Legacy DOLFIN
solver = PETScKrylovSolver("gmres", "ilu")
# or via options database
PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "ilu")
```
```python
# 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:
```python
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. DIVERGED_PC_FAILED: 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 `Function`s 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".