Semi-analytic Jacobian for 1D flames by speth · Pull Request #2139 · Cantera/cantera

speth · 2026-06-16T17:15:46Z

Changes proposed in this pull request

The 1D flame solver currently computes a finite difference Jacobian for use with the Newton's method solver. The main efficiency this exploits is the block tridiagonal pattern of the Jacobian, where the state at each grid point only affects the residuals at that grid point and the immediate neighbors, so everything scales linearly with the number of grid points $N$. However, the number of function evaluations to compute each block row is $O(K)$ and the cost is also $O(K)$, making this Jacobian evaluation an $O(NK^2)$ cost. Using a banded matrix solver, where the bandwidth is proportional to $K$, the factorization cost is $O(NK^3)$ and the backsubstitution/solve cost is $O(NK^2)$. It turns out, though, that the constant factors on these different operations are quite different, and for a range of mechanism sizes, the Jacobian construction cost dominates.

This cost can be dramatically reduced by computing a semi-analytical Jacobian, where the derivatives with respect to species mass fractions are computed analytically and derivatives with respect to the remaining state variables are computed using the existing finite difference approach. Derivatives of the reaction rate expressions make use of the existing Kinetics derivatives, while derivatives of the transport-related terms are new additions to the Flow1D class.

Currently, the analytic Jacobian works for ideal gases using mixture-averaged transport, including the mixture-averaged approximation for the Soret diffusion term. It is enabled by default, but automatically falls back to the finite difference form for other flow domains (specifically, IonFlow), transport models, or equations of state. It can also be manually switched via the Domain1D.jacobian_mode property.

This PR also revises the performance statistics available from the 1D solver, introducing an AnyMap-based store accessible via Sim1D.solver_stats that provides break downs of run time and call count for function evaluations, Jacobian evaluations, Jacobian factorizations, and linear solves for each grid adaptation step.

If applicable, provide an example illustrating new features this pull request is introducing

The following is a performance comparison on a laminar flame speed calculation, using mechanisms of various sizes, including a couple from NUI Galway converted to YAML format.

import cantera as ct
import pandas as pd
ct.suppress_thermo_warnings()

def run(mech, reactants, jacobian_mode):
      p = ct.one_atm  # pressure [Pa]
      Tin = 300.0  # unburned gas temperature [K]
      width = 0.03  # m
      loglevel = 0  # amount of diagnostic output (0 to 8)

      gas = ct.Solution(mech)
      gas.TPX = Tin, p, reactants

      f = ct.FreeFlame(gas, width=width)
      f.flame.jacobian_mode = jacobian_mode
      f.set_refine_criteria(ratio=3, slope=0.06, curve=0.12)
      f.transport_model = 'mixture-averaged'
      f.flux_gradient_basis = "mass"
      f.solve(loglevel=loglevel, auto=True)
      print(f"mixture-averaged flamespeed = {f.velocity[0]:7f} m/s")

      f.soret_enabled = True
      f.solve(loglevel=loglevel) # don't use 'auto' on subsequent solves
      print("mixture-averaged flamespeed with Soret diffusion = {f.velocity[0]:7f} m/s")

      df = pd.DataFrame.from_dict(f.solver_stats)
      print(df)
      print(df[['residual_time', 'jacobian_time', 'factor_time', 'solve_time', 'total_time']].sum())


CASES = [
    ('h2o2.yaml', 'H2:1.1, O2:1, AR:5'),  # 10 species, 29 reactions
    ('gri30.yaml', 'CH4:0.45, O2:1, N2:3.76'),  # 53 species, 325 reactions
    ('methanol-galway2016.yaml', 'CH3OH:0.666, O2:1, N2:3.76'),  # 173 species, 1011 reactions
    ('ethanol-DME-galway2018.yaml', 'C2H5OH:0.3666, O2:1, N2:3.76')  # 433 species, 1004 reactions
]

for mech, reactants in CASES:
      for jacobian_mode in ['finite-difference', 'analytic']:
            print("-" * 80)
            print(f"{mech=}")
            run(mech, reactants, jacobian_mode)
            print(f"{jacobian_mode=}")

Collating the output shows the resulting performance benefits (times in seconds; using OpenBLAS with OMP_NUM_THREADS=1).

n_species	mode	f eval	J eval	J factor	lin solve	total	speedup
10	FD	0.07	0.19	0.02	0.03	0.32	-
10	analytic	0.07	0.09	0.02	0.03	0.22	1.44
53	FD	1.21	9.62	0.55	0.34	11.76	-
53	analytic	1.17	1.89	0.54	0.32	3.95	2.97
173	FD	20.79	218.42	18.05	12.27	270.30	-
173	analytic	18.60	21.95	16.32	10.91	68.48	3.95
433	FD	92.42	632.96	200.89	62.40	992.94	-
433	analytic	96.15	69.33	196.78	65.21	431.70	2.30

The benefit for very small mechanisms is limited, but grows for medium-sized mechanisms. Eventually, for larger mechanisms, even though the Jacobian evaluation cost is reduced by a factor of ~9, the factorization step becomes the dominant contributor to the overall runtime and the relative speedup starts to decline.

AI Statement (required)

Extensive use of generative AI. Significant portions of code or documentation were generated with AI, including
logic and implementation decisions. All generated code and documentation were reviewed and understood by the contributor. Implemented using Claude Code (Opus 4.8) and superpowers.

Checklist

The pull request includes a clear description of this code change
Commit messages have short titles and reference relevant issues
Build passes (scons build & scons test) and unit tests address code coverage
Style & formatting of contributed code follows contributing guidelines
AI Statement is included
The pull request is ready for review

…o Python Add `Sim1D.eval_jacobian()` to trigger `evalSSJacobian()` from Python, and `Domain1D.global_component_index(name, j)` to map a component/point pair to its index in the global solution vector. Both are needed by the upcoming analytic Jacobian comparison harness. Includes a test that verifies the assembled Jacobian matches a manual finite-difference residual perturbation to within 2% (the tolerance accounts for frozen transport properties in the solver FD vs full recomputation in sim.eval()). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Adds a per-domain `jacobian_mode` flag ("finite-difference" or "analytic"), two virtual hooks (`hasAnalyticJacobian(j,n)` and `evalJacobianAnalytic(x, jac)`), and the column-skipping logic in `OneDim::evalJacobian` that calls `evalJacobianAnalytic` after the FD loop. No domain claims any columns yet, so this is a behavioral no-op. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Add an overload netProductionRates_ddCi(SparseMatrix&) that fills a caller-owned sparse matrix and reuses its nonzero pattern across calls. The pattern of d(wdot_k)/d(c_m) is fixed for a given mechanism and derivative settings, so it is built once (delegating to the allocating path) and subsequent calls overwrite the value array in place with no allocation or symbolic SpGEMM. The BulkKinetics override fuses the stoichiometry-matrix product directly into the compressed column storage using an O(K) dense scatter column (m_rbuf3) rather than a K x K dense buffer, mirroring calculateCompositionDerivatives() with ddX=false. Result is numerically identical to the allocating netProductionRates_ddCi() (verified to <1e-12 relative). No new Python API. Benchmark (gri30, optimize=n library): ~1447 us/call allocating vs ~883 us/call reused = 1.64x on the primitive. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>