Fix MultiPhaseEquil constant HP/SP solve interaction with temperature bounds

src/equil/MultiPhase.cpp

@@ -458,13 +458,16 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, double err, int maxsteps,
         double Thigh = 2.0*m_Tmax; // upper bound on T
         double Hlow = Undef, Hhigh = Undef;
         for (int n = 0; n < maxiter; n++) {
-            // if 'strt' is false, the current composition will be used as
-            // the starting estimate; otherwise it will be estimated
-            MultiPhaseEquil e(this, strt);
-            // start with a loose error tolerance, but tighten it as we get
-            // close to the final temperature
-
             try {
+                // if 'strt' is false, the current composition will be used as the
+                // starting estimate; otherwise it will be estimated. Construction can
+                // fail if the temperature guess puts a phase outside its valid
+                // temperature range while it still has moles in the current
+                // composition; the catch block below recovers by recomputing the
+                // initial composition (strt=true).
+                MultiPhaseEquil e(this, strt);
+                // start with a loose error tolerance, but tighten it as we get close to
+                // the final temperature
                 e.equilibrate(TP, err, maxsteps, loglevel);
                 double hnow = enthalpy();
                 // the equilibrium enthalpy monotonically increases with T;
@@ -533,9 +536,12 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, double err, int maxsteps,
         double Tlow = 1.0; // lower bound on T
         double Thigh = 1.0e6; // upper bound on T
         for (int n = 0; n < maxiter; n++) {
-            MultiPhaseEquil e(this, strt);
-
             try {
+                // Construction (inside the try so the catch below can recover) can fail
+                // if the temperature guess puts a condensed phase outside its valid
+                // range while it still has moles in the current composition; recovery
+                // recomputes the initial composition (strt=true).
+                MultiPhaseEquil e(this, strt);
                 e.equilibrate(TP, err, maxsteps, loglevel);
                 double snow = entropy();
                 if (snow < s0) {

src/equil/MultiPhaseEquil.cpp

@@ -7,6 +7,7 @@
 
 #include "cantera/equil/MultiPhaseEquil.h"
 #include "cantera/thermo/MolalityVPSSTP.h"
+#include "cantera/numerics/eigen_dense.h"
 #include "cantera/base/stringUtils.h"
 
 #include <cstdio>
@@ -72,11 +73,10 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
     //
     // When start=true (the default for TP equilibration), the solver computes
     // its own initial composition from elemental totals. In that case, if an
-    // excluded species has non-zero initial moles, its elemental contribution
-    // is tracked in b_missing and redistributed to valid component species
-    // before calling setInitialMoles.
-    vector<double> b_missing(m_nel, 0.0);
-    bool has_excluded_moles = false;
+    // excluded species has non-zero initial moles, it is recorded here so that
+    // its elemental contribution can be redistributed to valid component
+    // species before calling setInitialMoles.
+    vector<size_t> excluded_with_nonzero_moles;
     for (size_t k = 0; k < m_nsp_mix; k++) {
         size_t ip = m_mix->speciesPhaseIndex(k);
         if (!m_mix->solutionSpecies(k) &&
@@ -85,15 +85,11 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
             if (m_mix->speciesMoles(k) > 0.0) {
                 if (!start) {
                     throw CanteraError("MultiPhaseEquil::MultiPhaseEquil",
-                        "condensed-phase species {} is excluded since its thermo "
-                        "properties are\nnot valid at this temperature, but it has "
-                        "non-zero moles in the initial state.", m_mix->speciesName(k));
+                        "Species {} is excluded since its thermo properties are not "
+                        "valid\nat this temperature, but it has non-zero moles in the "
+                        "initial state.", m_mix->speciesName(k));
                 }
-                for (size_t m = 0; m < m_nel; m++) {
-                    b_missing[m] += m_mix->speciesMoles(k) *
-                                    m_mix->nAtoms(k, m_element[m]);
-                }
-                has_excluded_moles = true;
+                excluded_with_nonzero_moles.push_back(k);
             }
         }
     }
@@ -137,14 +133,32 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
     // linear Gibbs minimization. In this case, only the elemental composition
     // of the initial mixture state matters.
     if (start) {
-        if (has_excluded_moles) {
-            // Adjust m_moles to restore element contributions from excluded
-            // condensed-phase species. After Gaussian elimination, m_A has an identity
-            // matrix in the component-species columns, so adding b_missing[m] to
-            // m_moles[m_order[m]] changes only element m's total.
+        if (!excluded_with_nonzero_moles.empty()) {
+            // Adjust m_moles to restore the element contributions of excluded
+            // condensed-phase species. computeN() selects a set of linearly independent
+            // component species; the change in element totals from adding moles of
+            // these components is B*delta, where B is the component formula matrix
+            // B(i,j) = nAtoms(component_j, element_i). To restore exactly the missing
+            // element totals b_missing, we solve B*delta = b_missing and add delta[j]
+            // moles of component j. Note that b_missing must be computed after
+            // computeN(), since computeN() may reorder m_element or reduce m_nel.
             computeN();
-            for (size_t m = 0; m < m_nel; m++) {
-                m_moles[m_order[m]] += b_missing[m];
+            Eigen::VectorXd b_missing = Eigen::VectorXd::Zero(m_nel);
+            for (size_t k : excluded_with_nonzero_moles) {
+                for (size_t m = 0; m < m_nel; m++) {
+                    b_missing[m] += m_mix->speciesMoles(k) *
+                                    m_mix->nAtoms(k, m_element[m]);
+                }
+            }
+            Eigen::MatrixXd B(m_nel, m_nel);
+            for (size_t i = 0; i < m_nel; i++) {
+                for (size_t j = 0; j < m_nel; j++) {
+                    B(i, j) = m_mix->nAtoms(m_species[m_order[j]], m_element[i]);
+                }
+            }
+            Eigen::VectorXd delta = B.colPivHouseholderQr().solve(b_missing);
+            for (size_t j = 0; j < m_nel; j++) {
+                m_moles[m_order[j]] += delta[j];
             }
             updateMixMoles();
         }

test/python/test_equilibrium.py

@@ -619,3 +619,61 @@ def test_excluded_species_extra_element(solver, include_h2o):
     assert mix.phase_moles("PbO_yw") == approx(0.0)
     assert mix.phase_moles("PbO2_s") == approx(1.0, rel=1e-6)
     assert mix.phase_moles() == approx(phase_moles_ref, rel=1e-6)
+
+
+def _methane_air_water_mixture(T):
+    """
+    Stoichiometric methane/air mixture with excess liquid water as a heat sink, at
+    temperature ``T`` and 1 atm. Used by the excluded-species regression tests below.
+    """
+    liquid = ct.Solution('water.yaml', 'liquid_water')
+    gas = ct.Solution('gri30.yaml', transport_model=None)
+    gas.TPX = 300, ct.one_atm, 'O2:1.0, N2:3.76, CH4:0.5'
+    mix = ct.Mixture([(gas, 1.0), (liquid, 2.0)])
+    mix.T = T
+    mix.P = ct.one_atm
+    return mix
+
+
+@pytest.mark.parametrize("solver", ["vcs", "gibbs"])
+def test_excluded_species_HP(solver):
+    """
+    Constant enthalpy/pressure equilibration of a methane/air mixture with excess liquid
+    water acting as a heat sink. The temperature search visits temperatures where liquid
+    water's thermo data are invalid; the solver must exclude liquid water at those
+    temperatures while redistributing its (multi-element) elemental content to the
+    polyatomic component species without corrupting the element abundances. Regression
+    test for https://github.com/Cantera/cantera/issues/268.
+    """
+    mix = _methane_air_water_mixture(300)
+    elements = ['H', 'O', 'C', 'N']
+    elements_before = [mix.element_moles(e) for e in elements]
+
+    mix.equilibrate('HP', solver=solver)
+
+    # Element abundances must be conserved through the temperature search
+    assert [mix.element_moles(e) for e in elements] == approx(elements_before, rel=1e-7)
+    assert_at_equilibrium(mix)
+
+
+@pytest.mark.parametrize("solver", ["vcs", "gibbs"])
+def test_excluded_species_TP_polyatomic(solver):
+    """
+    Constant temperature/pressure equilibration at a temperature where liquid water's
+    thermo data are invalid, but with liquid water present in the initial composition.
+    The solver excludes liquid water and must redistribute its multi-element (H, O)
+    content onto the *polyatomic* gas-phase component species. This is the narrowest
+    direct reproducer of the redistribution bug: when the component formula matrix is not
+    the identity (i.e. components are not monatomic), adding the missing element totals
+    directly to the component moles corrupts the element abundances. Regression test for
+    the redistribution introduced in https://github.com/Cantera/cantera/pull/2116.
+    """
+    mix = _methane_air_water_mixture(1000.0)  # liquid water is excluded at this T
+    elements = ['H', 'O', 'C', 'N']
+    elements_before = [mix.element_moles(e) for e in elements]
+
+    mix.equilibrate('TP', solver=solver)
+
+    # Element abundances must be conserved despite the excluded liquid water
+    assert [mix.element_moles(e) for e in elements] == approx(elements_before, rel=1e-7)
+    assert_at_equilibrium(mix)