Add linear, quadratic, and geometric energy-grid options for the Boltzmann
two-term EEDF solver, selectable from the phase YAML `electron-energy-distribution`
node, in addition to the already-present custom grid.
The `Boltzmann-two-term` distribution type now honors a user-supplied
`energy-levels` grid via setCustomGrid instead of silently falling back to the
hard-coded default grid.
Relax test_eedf_solver to verify that the custom grid is honored and that the
resulting EEDF is a finite, positive, monotonically decreasing, normalized
distribution.
Kinetics derivatives computed using internal finite differencing
(for example, ddM for third-body reactions and pressure derivatives for
P-log and Chebyshev reactions) used an inconsistent multiplier state. This
fixes these calculations to use compute the difference based on the
unperturbed forward rate constants.
IonFlow overrides the species residual and diffusive-flux physics (ion drift,
electric-field coupling) without supplying matching analytic Jacobian
derivatives, so the inherited Flow1D analytic columns describe a different
residual than IonFlow solves. With the analytic Jacobian now the default, this
broke the stage-two ion-flame solve. Add an analyticJacobianSupported() hook
(false for IonFlow) so "auto" silently uses finite differences and explicit
"analytic" raises for ion flames.
Also pin finite differences in the high-pressure regrid test helper, whose
exact step-count assertions depend on the Jacobian method.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
In explicit jacobian_mode = "analytic", OneDim::evalJacobian now validates each
domain before the finite-difference column loop and raises if the kinetics lacks
composition derivatives or multicomponent transport is active. Internal/transient
conditions (adjoint force-full-update, fewer than 3 points) still fall back
silently.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Introduces a three-valued jacobian_mode. 'auto' (the new default) uses the
analytic Jacobian where supported and silently falls back to finite
differences otherwise; the capability probe no longer warns on fallback.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
accumulateStoichProduct now scans for product entries that the
supplied pattern cannot store and throws instead of silently
dropping them (and leaking into the next column). Compiled out
under NDEBUG.
Aso, merge the three analytic-vs-FD comparison classes into one
TestAnalyticVsFD with compare_modes as a member; solve the comparison
configs with auto=False to avoid masking unexpected failures; fold the
deprecated-alias check into test_solver_stats; and replace GRI 3.0 with
h2o2 (loose refinement) in test_full_solve_matches_fd.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Key changes:
- Extend analytic Jacobian to cover grid points 1 and N-2 (previously 2..N-3);
boundary residuals have no dependence on neighboring species so those rows are
safely skipped. Minimum grid size drops from 6 to 3.
- Refactor fluxJacobian to template<bool AtQ> computing only the needed endpoint
derivative, eliminating the scratch buffer and halving the inner-loop work.
- Remove unused x parameter from computeWdotDerivatives.
- Rename fuseStoichProduct → accumulateStoichProduct; explain template necessity
(StoichManagerN::derivatives returns Eigen::Map, not SparseMatrix).
- Move factorizeJacobian() body to SteadyStateSystem.cpp, localizing <chrono>.
- Inline the single-use accumulateGridTime lambda in SteadyStateSystem::solve.
- Reorganize Flow1D::hasAnalyticJacobian/evalJacobianAnalytic into the main
public section rather than a mid-class public: block.
- Fix math-mode formula in radiationPolyFactor docstring.
- Fix Doxygen group comment in SteadyStateSystem.h staging accumulators.
- Remove unused Cython pxd getters (gridSizeStats etc.); add pyi stubs for
new API (jacobian_mode, global_component_index, eval_jacobian, solver_stats).
- Add analytic-jacobian.md derivation page to Sphinx docs.
- Update release note with mechanism-size scaling explanation.
- Speed up test_radiation by using h2o2.yaml instead of GRI 3.0.
- Update TestAnalyticVsFD claimed-column range to match new j >= 1 condition.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Add TestAnalyticSolve.test_full_solve_matches_fd which solves a GRI-3.0
stoichiometric CH4/air free flame twice — once with the default
finite-difference Jacobian and once with analytic mode — exercising
regridding, time stepping, and workspace reallocation. Asserts flame
speeds agree to within 1e-4 relative tolerance (actual: ~7.6e-7).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add TestAnalyticConfigMatrix with four tests covering: axisymmetric
counterflow diffusion flame, unstrained burner flame, Soret diffusion
enabled, and two-point control. All pass without new C++ changes,
confirming Tasks 1-7 cover these row/column variants correctly.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds the compare_modes harness and TestAnalyticVsFD, and fixes two bugs
found by the FD oracle:
- The capability probe ran inside evalJacobianAnalytic (after the FD
column-skip loop), so claimed columns were never skipped and the
analytic post-pass double-counted them. The probe now runs lazily from
the const hasAnalyticJacobian()/usingAnalyticJacobian() queries.
- The 1/(rho cp) chain term in the energy row had the wrong sign
(FD-verified: d(1/(rho cp))/dY_m = (1/(rho cp))(wbar/W_m - cpm/cp),
giving -B*rcp*(...) for the residual rsd = -u dTdz - B/(rho cp)).
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>
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 `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>
#2134 introduced new thermodynamic consistency tests.
Skip them are there are not defined for two-temperature systems, or because the relative error is a bit too big for the acceptable error.
Add tests checking for consistency with the following identities:
- dG composition <-> T Maxwell relation: s_k = - d(mu_k)/dT at constant P and X
- dG composition <=> P Maxwell relation: v_k = d(mu_k)/dP at constant T and X
- dG pressure coefficiient: V = dG/dP at constant T
- dH pressure coefficient: dH/dP at constant T = V - T*dV/dT at constant P
- composition Hessian symmetry: d(mu_k)/dN_j = d(mu_j)/dN_k
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add a `gibbs_duhem_const_T_P` test to the parametrized ThermoPhase consistency
suite. It perturbs composition at constant T and P and checks the differential
identity sum_k X_k * dmu_k = 0, the complement to the Euler relation
g = sum_k X_k * mu_k that the suite already checked. This is the identity that
was missing: all existing tests held composition fixed and only perturbed T and
P, so a model whose chemical potentials are not consistent partial molar Gibbs
energies of a single G(T, P, N) could pass every check.
The finite-difference step is scaled to the smaller of the two perturbed mole
fractions so that the truncation error is independent of how dilute the species
are; a fixed absolute step gives false positives on dilute aqueous solutes.
Mark the Debye-Huckel B-dot model with per-ion size parameters and the
coverage-dependent surface model as known failures, since both report chemical
potentials that violate Gibbs-Duhem.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Fixes issue #2067 while also making the method const. After this
update, the state of the associated Solution is no longer updated,
which is a clear separation from the alternative setState method.
The previous handling of this redistribution (introduced in #2116 / 2219f65)
did not correctly handle polyatomic species.
Partially addresses #268
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The optimum-basis check at the end of each VCS iteration can disagree with
the basis that vcs_basopt actually selects: the check ranks a noncomponent
against a component by mole number alone, while vcs_basopt additionally
enforces linear independence and uses slightly different size weighting.
When they disagree, the check forces a component recalculation that re-selects
the same basis, which immediately trips the check again. The solver cycles
through its main loop forcing a new basis every iteration and never reaches
the equilibrium check, so it cannot converge.
Fixed by counting consecutive recalculations forced by this check and stopping
forcing once they clearly exceed what a real basis improvement requires. The
convergence criterion is the same for any equivalent basis, so suppressing
an oscillating swap is safe: it affects only conditioning, not the result.
Add a regression test covering a representative spread of the affected
compositions.
Fixes#266.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Previously, computeN only reselected the basis when the most-abundant species
per element wasn't a component. It never noticed when a component itself became
depleted. So a radical sitting in the basis at ~8e-10 mol could stay in the basis
and throttle every reaction that must consume it, since the step is bounded by
moles/|Δmoles|. The solver then settles into a limit cycle and never converges,
This is fixed by forcing component reselection when any component drops below a
threshold of 1e-3 * total_moles. This is safe because getComponents is idempotent-optimal;
it re-selects the same species if it's genuinely the best linearly-independent choice.
Otherwise, it swaps in a more abundant species.
Fixes#265
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Ray Speth
Tim E. Dawson
Ingmar Schoegl
Gaetanosaure
pag1pag
Christopher Neal
