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

Diagnose linear-solver failures in FEniCS/DOLFINx from the PETSc converged reason: DIVERGED_ITS (tolerances too tight, CG on non-SPD matrix, singular system needing a nullspace, preconditioner scaling), DIVERGED_PC_FAILED (direct-solver memory, singular factorization), DIVERGED_NANORINF (NaN in assembly). Includes the pure-Neumann nullspace recipe and solver/preconditioner selection for legacy DOLFIN and DOLFINx.

Exact reference: {"kind":"skill_version","skill_id":"skl_XNhRURC2BGwgOIyir8klEw","version_id":"skv_URzUmE3DwH7VXkrTdG5UdA"}

Applicability: [{"constraint":"Any legacy dolfin release (e.g. 2019.1.0); PETScKrylovSolver parameters absolute_tolerance, relative_tolerance, maximum_iterations","technology":"Legacy FEniCS (DOLFIN) PETScKrylovSolver","version_scheme":"unknown"},{"constraint":">=0.6.0; applies when an iterative KSP is selected via petsc_options or built manually (LinearProblem defaults to direct LU)","technology":"FEniCSx (DOLFINx) PETSc KSP","version_scheme":"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`:

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


## Supporting basis and limitations

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

## Change and rationale

New skill: diagnose PETSc Krylov solver divergence in FEniCS from the PETSc reason code.

'Unable to solve linear system using PETSc Krylov solver' is a frequent FEniCS Discourse error whose single most informative token — the PETSc reason code — is routinely ignored, leading to blind solver/preconditioner swapping. The reason code separates three unrelated problems (iteration limits, preconditioner setup failure, NaN assembly), and the DIVERGED_ITS branch further separates tolerance, solver-matrix mismatch, and singular-system causes.