PyBaMM Solver Convergence Failures: Numerics or Physics

A long PyBaMM degradation run fails at cycle 300 with a wall of SUNDIALS output ending in something like Corrector convergence failed repeatedly or with |h| = hmin. The first instinct is to reach for tolerances or a smaller timestep. That instinct is usually wrong, and following it costs hours.

This article is about telling the two cases apart. A convergence failure is either a numerical problem or a physical one, and the fixes are opposite. Loosening tolerances on a physically impossible state does not rescue the run – it hides the point where the model stopped meaning anything and lets the solver produce numbers you will later trust.

1. What the error actually says

IDAKLU is a BDF solver: at each step it forms an implicit algebraic system and solves it with Newton iterations. The “corrector” is that Newton loop. Two things can end a step:

• The corrector does not converge. Newton iterated up to its limit without the residual falling far enough. The solver responds by shrinking the timestep and retrying.
• The timestep hits hmin. After repeated shrinking there is nothing left to shrink. This is what |h| = hmin means.

The critical implication is in that sequence: by the time you see this error, the solver has already tried making the step smaller, many times. Manually reducing the step size or the output interval repeats work the adaptive controller has done and failed at. If shrinking the step were going to help, it already would have.

So the message is not “the step was too big”. It is closer to “at this state, the Newton iteration cannot find a solution at any step size I am willing to take”.

2. The decision that matters: numerics or physics

Before touching a single solver option, establish which of these you have. The test is cheap.

import pybamm

# Re-run the same experiment on a simpler model. SPMe drops the spatial
# resolution of DFN but keeps electrolyte gradients, so it survives states
# that DFN cannot represent.
for cls in (pybamm.lithium_ion.SPMe, pybamm.lithium_ion.DFN):
    sim = pybamm.Simulation(cls(), parameter_values=params, experiment=experiment)
    try:
        sol = sim.solve()
        print(f"{cls.__name__}: completed {len(sol.cycles)} cycles")
    except Exception as exc:
        print(f"{cls.__name__}: failed - {type(exc).__name__}: {exc}")

Read the result like this:

In practice on long degradation runs the first row is the common one, and it is the one people spend longest fighting with the wrong tool.

3. Finding the physically impossible state

When the answer is physics, the question becomes which variable went out of bounds. The usual suspects are bounded quantities that the equations assume stay inside their range:

• Particle surface concentration hitting 0 or the maximum. Butler-Volmer kinetics contain terms in c_s_surf and (c_s_max - c_s_surf); drive either to zero and the exchange current density collapses or the overpotential diverges.
• Electrolyte concentration reaching zero somewhere in the separator or a pore, which makes conductivity and diffusivity undefined.
• Porosity driven to zero by SEI growth. This one is specific to long degradation runs, and it is why a simulation can complete 250 cycles and fail on the 251st.

The way to see it is to solve up to just before the failure and inspect, rather than to stare at the solver log:

import numpy as np

# Solve a shorter run that succeeds, then look at how close the bounded
# variables came to their limits. A run that "worked" but approached a
# boundary is the same failure a few dozen cycles earlier.
sol = sim.solve()

for name in ("Negative particle surface concentration",
             "Positive particle surface concentration",
             "Electrolyte concentration [mol.m-3]",
             "Negative electrode porosity"):
    try:
        v = sol[name].entries
    except KeyError:
        continue          # names vary with model options; skip rather than crash the diagnostic
    print(f"{name:<44} min={np.min(v):.4g}  max={np.max(v):.4g}")

A stoichiometry that reaches 0.999 or 0.001, an electrolyte concentration an order of magnitude below its initial value, or a porosity heading toward zero all point at the same conclusion: the experiment is asking for something the cell cannot do. The fix is in the protocol or the parameters, not the solver.

4. When it really is numerics: the options that exist

These are IDAKLU's documented defaults. Knowing the numbers matters, because most advice on forums adjusts them without saying what they were.

The trap: options that silently do nothing

The IDAKLU options dictionary only takes effect when model.convert_to_format == "casadi". Set them on a model in another format and they are accepted without complaint and quietly ignored. You will conclude the option "did not help" when it was never applied.

solver = pybamm.IDAKLUSolver(
    rtol=1e-6,
    atol=1e-8,
    options={"max_nonlinear_iterations": 100, "max_order_bdf": 3},
)

model = pybamm.lithium_ion.DFN()

# Assert rather than assume: if this is not "casadi", every option above is decoration
assert model.convert_to_format == "casadi", (
    f"IDAKLU options are ignored unless convert_to_format is 'casadi'; "
    f"got {model.convert_to_format!r}"
)

sim = pybamm.Simulation(model, parameter_values=params,
                        experiment=experiment, solver=solver)

5. Why long runs fail and short ones do not

A failure that appears only after hundreds of cycles is almost never a solver setting that was fine at cycle 3 and wrong at cycle 300. Something monotonic has been accumulating.

The mechanism is usually this: degradation submodels change the geometry the equations are solved on. SEI growth consumes porosity; lithium plating and active material loss shift the usable stoichiometry window. Each cycle the model is slightly closer to a boundary, and eventually a state that was merely stiff becomes singular.

Two consequences for how you debug:

• Reproduce the failure faster by accelerating degradation, not by running longer. Raise the SEI rate constant until the same failure appears at cycle 20, then diagnose there. A twenty-cycle reproduction is something you can iterate on.
• Record the cycle index in your failure reports. "Fails at cycle 300 of 500 with parameter set X" is diagnosable; "the simulation crashes" is not.

import pybamm

# Accelerate degradation so a cycle-300 failure reproduces inside 20 cycles.
# The goal is not a physically meaningful result - it is turning a multi-hour
# iteration loop into a multi-minute one.
fast = params.copy()
fast["SEI kinetic rate constant [m.s-1]"] *= 50

sim = pybamm.Simulation(
    pybamm.lithium_ion.DFN(options={"SEI": "solvent-diffusion limited"}),
    parameter_values=fast,
    experiment=pybamm.Experiment(
        [("Discharge at 1C until 3.0V", "Charge at 1C until 4.2V", "Hold at 4.2V until C/50")] * 30
    ),
)
try:
    sol = sim.solve()
    print(f"completed {len(sol.cycles)} cycles - raise the multiplier further")
except Exception as exc:
    print(f"reproduced: {type(exc).__name__}")

6. A checklist, in the order worth trying

1. Record where it failed. Cycle index, step within the cycle, and the last successful time.
2. Re-run on SPMe. If it also fails, stop tuning the solver - the problem is the parameters or the protocol.
3. Inspect bounded variables on the longest run that succeeds. Look for stoichiometry near 0 or 1, electrolyte concentration collapsing, porosity approaching zero.
4. Check the experiment against the cell. A cutoff voltage the degraded cell can no longer reach, or a C-rate it can no longer sustain, produces exactly this failure.
5. Only now, options. max_nonlinear_iterations first, then max_order_bdf, and assert convert_to_format == "casadi" before believing any of it took effect.
6. Re-verify the results after loosening anything. A run that completes because the acceptance threshold was relaxed is not automatically a run you can trust; compare a shorter segment against the tighter settings.

7. What to record when it works

A convergence failure that took a day to diagnose will happen again in six months on a different parameter set. What makes it cheap the second time is the record:

Model:        DFN, options {SEI: solvent-diffusion limited, ...}
Parameters:   OKane2022, modified: SEI kinetic rate constant = 1e-15
Experiment:   1C/1C CCCV, 500 cycles requested
Failure:      cycle 312, discharge step, t = ...
Diagnosis:    negative particle surface stoichiometry -> 0.9993
Root cause:   cutoff voltage unreachable after ~30% capacity fade
Fix:          cutoff moved to 2.8V; NOT a solver change
Solver:       IDAKLU, rtol 1e-4, atol 1e-6, defaults otherwise

Note the "NOT a solver change" line. Half the value of this record is stopping the next person - usually you - from reaching for tolerances again.

References

• PyBaMM IDAKLUSolver API reference
• PyBaMM discussion: corrector convergence failed repeatedly or with |h| = hmin
• PyBaMM discussion: understanding solver doesn't converge
• PyBaMM coupled degradation notebook

If the diagnosis lands on the parameters rather than the solver, continue with fitting PyBaMM parameters and identifiability.

See also: PyBaMM architecture and solver choice, on the SPM/SPMe/DFN trade-off and why IDAKLU suits long degradation runs.