From 9d461f3fe4fd33a92a7a09bfb6ede9174c6f480e Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 29 May 2026 08:52:22 -0400 Subject: [PATCH 1/3] [Equil] Fix trace-phase handling in MultiPhaseEquil::stepComposition Before computing the step, drop trace stoichiometric phases that are being consumed (zero their moles and their formation reaction's dxi, so element balance is preserved). Partially fixes #265 Co-Authored-By: Claude Opus 4.7 --- src/equil/MultiPhaseEquil.cpp | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index 5e489c1558..683d92cc0d 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -495,6 +495,31 @@ double MultiPhaseEquil::stepComposition(int loglevel) m_iter++; double grad0 = computeReactionSteps(m_dxi); + // Drop trace stoichiometric (single-species) condensed phases that the current step + // would consume. Such a phase, with a mole number near the numerical floor left + // over from initialization, would otherwise pin the step size omega to a + // vanishingly small value (omega <= -moles/deltaN in the loop below), stalling + // progress on every other reaction. Because the step is removing it (dxi < 0, i.e. + // dG/RT for its formation reaction is positive), the phase should be absent. We + // zero its mole number and suppress its formation reaction here. + double total_moles = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + total_moles += std::max(0.0, m_moles[k]); + } + double trace = Tiny * total_moles; + bool dropped = false; + for (size_t j = 0; j < nFree(); j++) { + size_t k = m_order[j + m_nel]; + if (!m_dsoln[k] && m_moles[k] > 0.0 && m_moles[k] < trace && m_dxi[j] < 0.0) { + m_dxi[j] = 0.0; + m_moles[k] = 0.0; + dropped = true; + } + } + if (dropped) { + updateMixMoles(); + } + // compute the mole fraction changes. if (nFree()) { multiply(m_N, m_dxi, m_work); From 76a119f9bb03b5c4a92664ec8a689d35402c8200 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 29 May 2026 21:36:20 -0400 Subject: [PATCH 2/3] [Equil] Use consistent default max_steps for equilibrate The Python interface only allowed 1000 steps by default, while the C++ interface allowed 50000. --- interfaces/cython/cantera/mixture.pyi | 2 +- interfaces/cython/cantera/mixture.pyx | 2 +- interfaces/cython/cantera/thermo.pyi | 2 +- interfaces/cython/cantera/thermo.pyx | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/interfaces/cython/cantera/mixture.pyi b/interfaces/cython/cantera/mixture.pyi index 260e2720a9..7f8f2714c5 100644 --- a/interfaces/cython/cantera/mixture.pyi +++ b/interfaces/cython/cantera/mixture.pyi @@ -61,7 +61,7 @@ class Mixture: XY: PropertyPair, solver: EquilibriumSolver = "auto", rtol: float = 1e-9, - max_steps: int = 1000, + max_steps: int = 50000, max_iter: int = 100, estimate_equil: Literal[-1, 0, 1] = 0, log_level: LogLevel = 0, diff --git a/interfaces/cython/cantera/mixture.pyx b/interfaces/cython/cantera/mixture.pyx index 84a7a77ad0..98bbd56e2e 100644 --- a/interfaces/cython/cantera/mixture.pyx +++ b/interfaces/cython/cantera/mixture.pyx @@ -286,7 +286,7 @@ cdef class Mixture: self.mix.getChemPotentials(view) return data - def equilibrate(self, XY, solver='auto', rtol=1e-9, max_steps=1000, + def equilibrate(self, XY, solver='auto', rtol=1e-9, max_steps=50000, max_iter=100, estimate_equil=0, log_level=0): """ Set to a state of chemical equilibrium holding property pair ``XY`` diff --git a/interfaces/cython/cantera/thermo.pyi b/interfaces/cython/cantera/thermo.pyi index 799cc4155e..952611c68c 100644 --- a/interfaces/cython/cantera/thermo.pyi +++ b/interfaces/cython/cantera/thermo.pyi @@ -162,7 +162,7 @@ class ThermoPhase(_SolutionBase): XY: PropertyPair, solver: EquilibriumSolver = "auto", rtol: float = 1e-9, - max_steps: int = 1000, + max_steps: int = 50000, max_iter: int = 100, estimate_equil: int = 0, log_level: LogLevel = 0, diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index d504e44f30..b13a608f54 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -408,7 +408,7 @@ cdef class ThermoPhase(_SolutionBase): return 1.0 def equilibrate(self, XY, solver='auto', double rtol=1e-9, - int max_steps=1000, int max_iter=100, int estimate_equil=0, + int max_steps=50000, int max_iter=100, int estimate_equil=0, int log_level=0): """ Set to a state of chemical equilibrium holding property pair From a9c0eff1015445858c74151e28465b74bd7f8cd7 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 29 May 2026 21:45:29 -0400 Subject: [PATCH 3/3] [Equil] Fix component basis reselection in MultiPhaseEquil::computeN MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- src/equil/MultiPhaseEquil.cpp | 33 ++++++-- test/python/test_equilibrium.py | 144 +++++++++++++++++++++++--------- 2 files changed, 131 insertions(+), 46 deletions(-) diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index 683d92cc0d..c6e7e8b872 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -709,7 +709,8 @@ void MultiPhaseEquil::computeN() m_sortindex[k] = moleFractions[k].second; } - for (size_t m = 0; m < m_nel; m++) { + bool reselect = m_force; + for (size_t m = 0; m < m_nel && !reselect; m++) { size_t k = 0; for (size_t ik = 0; ik < m_nsp; ik++) { k = m_sortindex[ik]; @@ -723,11 +724,33 @@ void MultiPhaseEquil::computeN() ok = true; } } - if (!ok || m_force) { - getComponents(m_sortindex); - m_force = true; - break; + if (!ok) { + reselect = true; + } + } + + // Reselect the component basis if any current component has been depleted to a + // negligible mole number. A near-empty basis species throttles every reaction that + // must consume it (the step size is bounded by moles/|deltaN|), which stalls the + // solver in a limit cycle without converging. Reselecting from the + // mole-fraction-sorted list replaces the depleted species with an abundant, + // linearly-independent one when possible. + if (!reselect) { + double total_moles = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + total_moles += std::max(0.0, m_moles[k]); } + for (size_t m = 0; m < m_nel; m++) { + if (m_moles[m_order[m]] < 1.0e-3 * total_moles) { + reselect = true; + break; + } + } + } + + if (reselect) { + getComponents(m_sortindex); + m_force = true; } } diff --git a/test/python/test_equilibrium.py b/test/python/test_equilibrium.py index 1c026272a5..22103682c6 100644 --- a/test/python/test_equilibrium.py +++ b/test/python/test_equilibrium.py @@ -9,6 +9,47 @@ compare ) + +def assert_at_equilibrium(mix, abs_tol=1e-6): + """ + Verify that ``mix`` is at chemical equilibrium by checking that the chemical + potential of each species satisfies ``mu_k = sum_e n_{e,k} * lambda_e`` where + ``lambda_e`` is the element potential for element e and ``n_{e,k}`` is the + number of atoms of element e in species k. We solve for ``lambda`` by weighted + least squares (weighting by moles so trace species do not dominate the fit), + then require ``|mu_k - sum_e n_{e,k} lambda_e| / (R T) < abs_tol`` for every + species whose moles are more than 1e-10 of the total. Returns the worst relative + residual. + """ + R = ct.gas_constant + T = mix.T + n_species = mix.n_species + n_elements = mix.n_elements + mu = np.asarray(mix.chemical_potentials) + moles = np.asarray(mix.species_moles) + composition = np.array([[mix.n_atoms(k, e) for e in range(n_elements)] + for k in range(n_species)]) + + w = np.sqrt(np.maximum(moles, 0.0)) + mask = w > 0.0 + Aw = composition[mask] * w[mask, None] + muw = mu[mask] * w[mask] + lam, *_ = np.linalg.lstsq(Aw, muw, rcond=None) + + total = moles.sum() + worst = 0.0 + for k in range(n_species): + if moles[k] > 1e-10 * total: + residual = (mu[k] - composition[k] @ lam) / (R * T) + assert abs(residual) < abs_tol, ( + f"species {mix.species_name(k)} (k={k}): " + f"moles={moles[k]:.3g} of total {total:.3g}, " + f"|mu_k - sum_e n_{{e,k}} lambda_e| / RT = {abs(residual):.3g}" + ) + worst = max(worst, abs(residual)) + return worst + + class EquilTestCases: """ Base class for equilibrium test cases, parameterized by the solver to use. @@ -180,46 +221,6 @@ def _h2o2_phase2_samples(n=20): stride = max(1, len(combos) // n) return [list(c) for c in combos[::stride][:n]] - @staticmethod - def _at_equilibrium(mix, abs_tol=1e-6): - """ - Verify that ``mix`` is at chemical equilibrium by checking that the chemical - potential of each species satisfies ``mu_k = sum_e n_{e,k} * lambda_e`` where - ``lambda_e`` is the element potential for element e and ``n_{e,k}`` is the - number of atoms of element e in species k. We solve for ``lambda`` by weighted - least squares (weighting by moles so trace species do not dominate the fit), - then require ``|mu_k - sum_e n_{e,k} lambda_e| / (R T) < abs_tol`` for every - species whose moles are more than 1e-10 of the total. Returns the worst relative - residual. - """ - R = ct.gas_constant - T = mix.T - n_species = mix.n_species - n_elements = mix.n_elements - mu = np.asarray(mix.chemical_potentials) - moles = np.asarray(mix.species_moles) - composition = np.array([[mix.n_atoms(k, e) for e in range(n_elements)] - for k in range(n_species)]) - - w = np.sqrt(np.maximum(moles, 0.0)) - mask = w > 0.0 - Aw = composition[mask] * w[mask, None] - muw = mu[mask] * w[mask] - lam, *_ = np.linalg.lstsq(Aw, muw, rcond=None) - - total = moles.sum() - worst = 0.0 - for k in range(n_species): - if moles[k] > 1e-10 * total: - residual = (mu[k] - composition[k] @ lam) / (R * T) - assert abs(residual) < abs_tol, ( - f"species {mix.species_name(k)} (k={k}): " - f"moles={moles[k]:.3g} of total {total:.3g}, " - f"|mu_k - sum_e n_{{e,k}} lambda_e| / RT = {abs(residual):.3g}" - ) - worst = max(worst, abs(residual)) - return worst - @pytest.mark.parametrize("phase2_names", _h2o2_phase2_samples()) @pytest.mark.parametrize("T,P_atm,N_inert", _H2O2_TPN) def test_equilibrate(self, T, P_atm, N_inert, phase2_names): @@ -247,7 +248,68 @@ def test_equilibrate(self, T, P_atm, N_inert, phase2_names): np.testing.assert_allclose(elements_after, elements_before, rtol=1e-7, atol=1e-12, err_msg="element moles changed during equilibration") - self._at_equilibrium(mix, abs_tol=1e-6) + assert_at_equilibrium(mix, abs_tol=1e-6) + + +class TestMultiphase_CHO: + """ + Regression tests for the 'gibbs' (``MultiPhaseEquil``) solver on C/H/O compositions + which previously failed to converge (Issue #265). Each case uses the GRI 3.0 gas + phase at 923 K and 1 atm, with the element amounts set by integer atom counts (C, H, + O). The cases exercise the two failure mechanisms that were fixed: + + - A trace single-species condensed phase (graphite) pinning the step size to a + negligible value (oxygen-rich cases such as C=23, H=79, O=98). + - A depleted species left in the component basis, which throttles every reaction + that must consume it and sends the solver into a non-converging limit cycle + (carbon-rich, low-oxygen cases such as C=77, H=120, O=3). + """ + + T = 923.0 + P = ct.one_atm + + @staticmethod + def _ids(cases): + return [f"C{C}-H{H}-O{O}" for C, H, O in cases] + + def _check(self, phases, C, H, O, max_steps): + gas = phases[0][0] + gas.TPX = self.T, self.P, {'C': C, 'H': H, 'O': O} + mix = ct.Mixture(phases) + elements_before = np.array([mix.element_moles(e) + for e in range(mix.n_elements)]) + mix.T = self.T + mix.P = self.P + mix.equilibrate('TP', solver='gibbs', max_steps=max_steps) + elements_after = np.array([mix.element_moles(e) + for e in range(mix.n_elements)]) + np.testing.assert_allclose( + elements_after, elements_before, rtol=1e-7, atol=1e-12, + err_msg="element moles changed during equilibration") + assert_at_equilibrium(mix, abs_tol=1e-6) + + # (C, H, O) atom counts. Single-phase (gas only) cases are carbon-rich, low-oxygen + # compositions that previously stalled in a component-basis limit cycle. + _single_phase = [ + (40, 147, 13), (77, 120, 3), (78, 119, 4), + (80, 117, 3), (77, 119, 4), (80, 116, 4), + ] + @pytest.mark.parametrize("C,H,O", _single_phase, ids=_ids(_single_phase)) + def test_single_phase(self, C, H, O): + gas = ct.Solution('gri30.yaml', transport_model=None) + self._check([(gas, 1.0)], C, H, O, max_steps=2000) + + # Two-phase (gas + graphite) cases: oxygen-rich (trace condensed-phase pinning) + # through carbon-rich (slow convergence with the condensed phase present). + _two_phase = [ + (23, 79, 98), (124, 67, 9), (160, 36, 4), + (150, 46, 4), (100, 90, 10), (40, 147, 13), + ] + @pytest.mark.parametrize("C,H,O", _two_phase, ids=_ids(_two_phase)) + def test_two_phase(self, C, H, O): + gas = ct.Solution('gri30.yaml', transport_model=None) + carbon = ct.Solution('graphite.yaml') + self._check([(gas, 1.0), (carbon, 0.0)], C, H, O, max_steps=5000) @pytest.fixture(scope='function')