Fix mole fractions near pure-species limit in BinarySolutionTabulated…

…Thermo Fixes behavior when the tabulated species mole fraction has become 1.0 to within machine precision, while the other species mole fraction is not quite zero.
Cantera · Mar 4, 2023 · c3eb05c · c3eb05c
1 parent f19bd3f
commit c3eb05c
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 13 deletions.
diff --git a/include/cantera/thermo/BinarySolutionTabulatedThermo.h b/include/cantera/thermo/BinarySolutionTabulatedThermo.h
@@ -240,9 +240,6 @@ class BinarySolutionTabulatedThermo : public IdealSolidSolnPhase
     //! Current tabulated species index
     size_t m_kk_tab;
 
-    //! Current tabulated species mole fraction
-    mutable double m_xlast;
-
     //! Tabulated contribution to h0[m_kk_tab] at the current composition
     mutable double m_h0_tab;
 

diff --git a/src/thermo/BinarySolutionTabulatedThermo.cpp b/src/thermo/BinarySolutionTabulatedThermo.cpp
@@ -22,8 +22,6 @@ namespace Cantera
 BinarySolutionTabulatedThermo::BinarySolutionTabulatedThermo(const std::string& inputFile,
                                                              const std::string& id_)
     : m_kk_tab(npos)
-    , m_xlast(-1)
-
 {
     initThermoFile(inputFile, id_);
 }
@@ -36,22 +34,24 @@ void BinarySolutionTabulatedThermo::compositionChanged()
 
 void BinarySolutionTabulatedThermo::_updateThermo() const
 {
-    double xnow = moleFraction(m_kk_tab);
-    bool x_changed = (m_xlast != xnow);
+    static const int cacheId = m_cache.getId();
+    CachedScalar cached = m_cache.getScalar(cacheId);
+    bool x_changed = !cached.validate(stateMFNumber());
 
     if (x_changed) {
-        m_h0_tab = interpolate(xnow, m_enthalpy_tab);
-        m_s0_tab = interpolate(xnow, m_entropy_tab);
-        if (xnow == 0) {
+        double x_tab = moleFraction(m_kk_tab);
+        double x_other = moleFraction(1 - m_kk_tab);
+        m_h0_tab = interpolate(x_tab, m_enthalpy_tab);
+        m_s0_tab = interpolate(x_tab, m_entropy_tab);
+        if (x_tab == 0) {
             m_s0_tab = -BigNumber;
-        } else if (xnow == 1) {
+        } else if (x_other == 0) {
             m_s0_tab = BigNumber;
         } else {
-            m_s0_tab += GasConstant*std::log(xnow/(1.0-xnow)) +
+            m_s0_tab += GasConstant*std::log(x_tab / x_other) +
                     GasConstant/Faraday*std::log(standardConcentration(1-m_kk_tab)
                     /standardConcentration(m_kk_tab));
         }
-        m_xlast = xnow;
     }
 
     double tnow = temperature();