From 25be1c8b0e344b9c15f043fba149e47afb6c5475 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 30 May 2026 15:31:52 -0400 Subject: [PATCH 1/3] [Equil] Enforce iteration limit in VCS solver main loop The VCS solver's main loop checked the iteration limit at the top of the MAIN stage but only set the failure flag (iconv = -1) without exiting the loop, relying on the subsequent equilibrium-check stage to act on the limit. That stage is not always reached: a component-basis swap oscillation can keep the solver cycling through the MAIN stage indefinitely, forcing a new component calculation on every iteration and never advancing to the equilibrium check. Fixed by routing the over-iteration case directly to the cleanup-and-return block. Co-Authored-By: Claude Opus 4.8 --- src/equil/vcs_solve_TP.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index f30d69423d5..281fe2d66b5 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -184,11 +184,15 @@ int VCS_SOLVE::solve_TP(int print_lvl, int printDetails, int maxit) forceComponentCalc = 0; iti = 0; } - // Check on too many iterations. If we have too many iterations, - // Clean up and exit code even though we haven't converged. - // -> we have run out of iterations! + // Check on too many iterations. If we have too many iterations, clean up + // and exit even though we haven't converged. The equilibrium check below + // also enforces the iteration limit, but it is not always reached: a + // basis-swap oscillation can keep the solver cycling through this MAIN + // stage indefinitely. if (m_VCount->Its > maxit) { iconv = -1; + stage = RETURN_A; + continue; } solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies, forceComponentCalc, stage, printDetails > 0); From e80536dda1b56688ff6004b79cb28078e31cc7ce Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 30 May 2026 15:32:37 -0400 Subject: [PATCH 2/3] [Equil] Filter VCS non-convergence messages based on print level Co-Authored-By: Claude Opus 4.8 --- src/equil/vcs_solve_TP.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 281fe2d66b5..39172c48267 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -311,11 +311,13 @@ int VCS_SOLVE::solve_TP(int print_lvl, int printDetails, int maxit) vcs_TCounters_report(); } - // FILL IN - if (iconv < 0) { - plogf("ERROR: FAILURE its = %d!\n", m_VCount->Its); - } else if (iconv == 1) { - plogf("WARNING: RANGE SPACE ERROR encountered\n"); + // Report non-convergence if diagnostic output has been requested. + if (print_lvl > 0) { + if (iconv < 0) { + plogf("ERROR: FAILURE its = %d!\n", m_VCount->Its); + } else if (iconv == 1) { + plogf("WARNING: RANGE SPACE ERROR encountered\n"); + } } // Return a Flag indicating whether convergence occurred return iconv; From a5ed78dc71396504fb5bf724c576d36e89968b15 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 30 May 2026 15:58:35 -0400 Subject: [PATCH 3/3] [Equil] Prevent component-basis oscillation in VCS solver 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 --- include/cantera/equil/vcs_solve.h | 2 +- src/equil/vcs_solve_TP.cpp | 91 +++++++++++++++++++++---------- test/python/test_equilibrium.py | 43 +++++++++++++++ 3 files changed, 105 insertions(+), 31 deletions(-) diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index bda899b57c1..64ddcca7ec2 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -784,7 +784,7 @@ class VCS_SOLVE void solve_tp_component_calc(bool& allMinorZeroedSpecies); void solve_tp_inner(size_t& iti, size_t& it1, bool& uptodate_minors, bool& allMinorZeroedSpecies, int& forceComponentCalc, - int& stage, bool printDetails); + int& stage, bool printDetails, int& basisOptCount); void solve_tp_equilib_check(bool& allMinorZeroedSpecies, bool& uptodate_minors, bool& giveUpOnElemAbund, int& solveFail, size_t& iti, size_t& it1, int maxit, diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 39172c48267..fce493e34f9 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -91,6 +91,12 @@ int VCS_SOLVE::solve_TP(int print_lvl, int printDetails, int maxit) bool uptodate_minors = true; int forceComponentCalc = 1; + // Counts consecutive component recalculations forced by the optimum-basis check in + // solve_tp_inner. A real basis improvement is resolved in a few recalculations; a + // large count means the check and vcs_basopt disagree on the best basis and are + // oscillating, so the forcing should be suppressed. + int basisOptCount = 0; + // Set the debug print lvl to the same as the print lvl. m_debug_print_lvl = printDetails; if (printDetails > 0 && print_lvl == 0) { @@ -194,8 +200,8 @@ int VCS_SOLVE::solve_TP(int print_lvl, int printDetails, int maxit) stage = RETURN_A; continue; } - solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies, - forceComponentCalc, stage, printDetails > 0); + solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies, + forceComponentCalc, stage, printDetails > 0, basisOptCount); lec = false; } else if (stage == EQUILIB_CHECK) { // EQUILIBRIUM CHECK FOR MAJOR SPECIES @@ -353,7 +359,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, bool& uptodate_minors, bool& allMinorZeroedSpecies, int& forceComponentCalc, - int& stage, bool printDetails) + int& stage, bool printDetails, + int& basisOptCount) { if (iti == 0) { // SET INITIAL VALUES FOR ITERATION @@ -1058,42 +1065,66 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } // CHECK FOR OPTIMUM BASIS - for (size_t i = 0; i < m_numRxnRdc; ++i) { - size_t k = m_indexRxnToSpecies[i]; - if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - continue; - } - for (size_t j = 0; j < m_numComponents; ++j) { - bool doSwap = false; - if (m_SSPhase[j]) { - doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) > - (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); - if (!m_SSPhase[k] && doSwap) { - doSwap = m_molNumSpecies_old[k] > (m_molNumSpecies_old[j] * 1.01); - } - } else { - if (m_SSPhase[k]) { + // + // This check can disagree with the basis actually selected by vcs_basopt: here a + // noncomponent is judged "better" than a component using mole numbers alone, while + // vcs_basopt additionally enforces linear independence and uses slightly different + // size weighting. When the two disagree, forcing a recalculation re-selects the + // same basis, which trips this check again, causing an oscillation that prevents + // the solver from ever reaching the equilibrium check. Because the convergence + // criterion does not depend on which equivalent basis is chosen, suppressing the + // swap once it is clearly oscillating is safe: it only affects conditioning, not + // the answer. Allow a run of recalculations to settle, then stop forcing. The + // threshold is set well above the number a legitimately converging case needs but + // far below the iteration limit, so a true oscillation is still broken promptly. + int maxBasisOpt = 200; + if (basisOptCount <= maxBasisOpt) { + for (size_t i = 0; i < m_numRxnRdc; ++i) { + size_t k = m_indexRxnToSpecies[i]; + if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + continue; + } + for (size_t j = 0; j < m_numComponents; ++j) { + bool doSwap = false; + if (m_SSPhase[j]) { doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); - if (!doSwap) { + if (!m_SSPhase[k] && doSwap) { doSwap = m_molNumSpecies_old[k] > (m_molNumSpecies_old[j] * 1.01); } } else { - doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) > - (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); + if (m_SSPhase[k]) { + doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) > + (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); + if (!doSwap) { + doSwap = m_molNumSpecies_old[k] > (m_molNumSpecies_old[j] * 1.01); + } + } else { + doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) > + (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); + } } - } - if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) { - if (m_debug_print_lvl >= 2) { - plogf(" --- Get a new basis because %s", m_speciesName[k]); - plogf(" is better than comp %s", m_speciesName[j]); - plogf(" and share nonzero stoic: %-9.1f\n", - m_stoichCoeffRxnMatrix(j,i)); + if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) { + if (m_debug_print_lvl >= 2) { + plogf(" --- Get a new basis because %s", m_speciesName[k]); + plogf(" is better than comp %s", m_speciesName[j]); + plogf(" and share nonzero stoic: %-9.1f\n", + m_stoichCoeffRxnMatrix(j,i)); + } + forceComponentCalc = 1; + basisOptCount++; + return; } - forceComponentCalc = 1; - return; } } + // The optimum-basis check found no improvement: the basis is genuinely settled, + // so reset the counter. When the check is suppressed because it is oscillating, + // the counter is deliberately left latched so the suppression persists for the + // rest of this solve. + basisOptCount = 0; + } else { + debuglog(" --- Optimum basis check suppressed: basis is oscillating\n", + m_debug_print_lvl >= 2); } debuglog(" --- Check for an optimum basis passed\n", m_debug_print_lvl >= 2); stage = EQUILIB_CHECK; diff --git a/test/python/test_equilibrium.py b/test/python/test_equilibrium.py index 22103682c6f..f22c0408ed1 100644 --- a/test/python/test_equilibrium.py +++ b/test/python/test_equilibrium.py @@ -312,6 +312,49 @@ def test_two_phase(self, C, H, O): self._check([(gas, 1.0), (carbon, 0.0)], C, H, O, max_steps=5000) +class TestMultiphase_VCS_CHO: + """ + Regression tests for the 'vcs' (``vcs_MultiPhaseEquil``) solver on C/H/O + compositions which previously failed to converge (Issue #266). Each case uses + the GRI 3.0 gas phase plus a graphite condensed phase at 923 K and 1 atm, with + the element amounts set by integer atom counts (C, H, O). + """ + + 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, C, H, O, max_steps): + gas = ct.Solution('gri30.yaml', transport_model=None) + carbon = ct.Solution('graphite.yaml') + gas.TPX = self.T, self.P, {'C': C, 'H': H, 'O': O} + mix = ct.Mixture([(gas, 1.0), (carbon, 0.0)]) + 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='vcs', 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) + + # A representative spread across the band of compositions (C=13..30) that + # previously triggered the basis-swap oscillation, including its edges. + _cases = [ + (13, 78, 9), (15, 70, 15), (16, 65, 19), (20, 60, 20), + (22, 56, 22), (25, 50, 25), (27, 44, 29), (30, 36, 34), + ] + @pytest.mark.parametrize("C,H,O", _cases, ids=_ids(_cases)) + def test_two_phase(self, C, H, O): + self._check(C, H, O, max_steps=10000) + + @pytest.fixture(scope='function') def extra_elements(request): s = """