Skip to content

Fix convergence failures in VCS constant TP solver #2124

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
speth wants to merge 3 commits into
base: main
Choose a base branch
from
+119 −39
Open

Fix convergence failures in VCS constant TP solver #2124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
25be1c8
[Equil] Enforce iteration limit in VCS solver main loop
speth May 30, 2026
e80536d
[Equil] Filter VCS non-convergence messages based on print level
speth May 30, 2026
a5ed78d
[Equil] Prevent component-basis oscillation in VCS solver
speth May 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion include/cantera/equil/vcs_solve.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
113 changes: 75 additions & 38 deletions src/equil/vcs_solve_TP.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
43 changes: 43 additions & 0 deletions test/python/test_equilibrium.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = """
Expand Down
Loading