diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index bda899b57c..64ddcca7ec 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 f30d69423d..fce493e34f 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) { @@ -184,14 +190,18 @@ 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); + 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 @@ -307,11 +317,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; @@ -347,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 @@ -1052,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 22103682c6..f22c0408ed 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 = """