Stiffness and Implicit Integration — Working Notes

Stiffness is one of those concepts that sounds abstract until you watch an explicit integrator explode on a hydrogen-air mechanism at 1200 K. Then it becomes very concrete.

What stiffness actually means

A system of ODEs is stiff when the ratio of the largest to smallest eigenvalue of the Jacobian is large — often many orders of magnitude in combustion. The practical consequence: explicit methods are forced to take timesteps governed by the fastest (most negative) eigenvalue even when the solution evolves on the timescale of the slowest one.

For a 9-species hydrogen mechanism, the fastest chemistry timescales are on the order of $10^{-10}$ s. The energy release and bulk flow evolve on $10^{-4}$ s or longer. An explicit solver must step at $\sim 10^{-10}$ s across an integration window of $10^{-3}$ s. That's $10^7$ steps for a problem that might need a few hundred implicit steps.

The implicit tradeoff

Implicit methods — BDF, SDIRK, and their kin — handle stiffness by solving a nonlinear system at each timestep instead of evaluating an explicit update. The timestep constraint switches from being eigenvalue-bound to accuracy-bound, which is the right regime.

The cost is that each timestep requires solving:

$F(y_{n+1}) = y_{n+1} - y_n - h \cdot f(t_{n+1}, y_{n+1}) = 0$

via Newton iterations. Each Newton iteration requires a Jacobian-vector product (or full Jacobian factorization for direct solvers). For large mechanisms ($N > 100$ species), this is where the real computational cost lives.

Jacobian structure in combustion

The species-energy Jacobian for a perfectly stirred reactor or 0-D ignition problem is dense in general, but has exploitable structure:

• Each species rate depends on a small subset of species (sparse reaction participation)
• Temperature sensitivity is often strongly coupled to a few radical species
• The Jacobian eigenvalue spectrum clusters: a few very stiff modes, many moderately stiff, a long tail of slow modes

Analytically computing the Jacobian (via cantera.Reactor.jacobian() or through pyJac) outperforms finite-difference approximations by 2–5× and gives cleaner structure for preconditioner construction.

Notes on Cantera usage

Cantera's ReactorNet uses CVODE (from SUNDIALS) under the hood, which is a BDF integrator with adaptive order and stepsize. A few things I keep coming back to:

• set_max_time_step() matters early in ignition: the solver sometimes takes overconfident steps before the radical pool builds up.
• set_initial_time() resets internal CVODE state — useful when chaining multiple reactor calculations without re-instantiating.
• For large mechanisms, switching from dense to banded or sparse Jacobian representation in CVODE cuts memory and time dramatically. Not exposed via the Python API by default; needs solver options.

Open questions I'm sitting with

• How much does Jacobian reuse (vs. recomputing each Newton step) matter for typical ignition problems? CVODE does this adaptively, but understanding the heuristic is worth digging into.
• Preconditioning strategy for GMRES in the context of Newton-Krylov: ILU is standard, but block-Jacobi with species-group clustering looks promising for larger mechanisms.
• Is there a clean way to expose the full Jacobian eigenvalue spectrum from a running CVODE integration? Would help with stiffness diagnostics.