Fix IdealSolidSolnPhase and BinarySolutionTabulatedThermo

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;
 
@@ -255,8 +252,6 @@ class BinarySolutionTabulatedThermo : public IdealSolidSolnPhase
     vector_fp m_entropy_tab;
     vector_fp m_molar_volume_tab;
     vector_fp m_derived_molar_volume_tab;
-    vector_fp m_partial_molar_volume_1_tab;
-    vector_fp m_partial_molar_volume_2_tab;
 
 private:
     virtual void _updateThermo() const;

include/cantera/thermo/IdealSolidSolnPhase.h

@@ -596,9 +596,10 @@ class IdealSolidSolnPhase : public ThermoPhase
 
     //! Vector of molar volumes for each species in the solution
     /**
-     * Species molar volumes \f$ m^3 kmol^-1 \f$
+     * Species molar volumes (\f$ m^3 kmol^-1 \f$) at the current mixture state.
+     * For the IdealSolidSolnPhase class, these are constant.
      */
-    vector_fp m_speciesMolarVolume;
+    mutable vector_fp m_speciesMolarVolume;
 
     //! Vector containing the species reference enthalpies at T = m_tlast
     mutable vector_fp m_h0_RT;

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();
@@ -127,8 +127,6 @@ void BinarySolutionTabulatedThermo::initThermo()
         m_entropy_tab.resize(N);
         m_molar_volume_tab.resize(N);
         m_derived_molar_volume_tab.resize(N);
-        m_partial_molar_volume_1_tab.resize(N);
-        m_partial_molar_volume_2_tab.resize(N);
 
         for (size_t i = 0; i < N; i++) {
             m_molefrac_tab[i] = x_h[i].first;
@@ -138,13 +136,6 @@ void BinarySolutionTabulatedThermo::initThermo()
         }
 
         diff(m_molar_volume_tab, m_derived_molar_volume_tab);
-
-        for (size_t i = 0; i < N; i++) {
-            m_partial_molar_volume_1_tab[i] = m_molar_volume_tab[i] +
-                    (1-m_molefrac_tab[i]) * m_derived_molar_volume_tab[i];
-            m_partial_molar_volume_2_tab[i] = m_molar_volume_tab[i] -
-                    m_molefrac_tab[i] * m_derived_molar_volume_tab[i];
-        }
     }
     IdealSolidSolnPhase::initThermo();
 }
@@ -162,6 +153,7 @@ void BinarySolutionTabulatedThermo::getParameters(AnyMap& phaseNode) const
     tabThermo["mole-fractions"] = m_molefrac_tab;
     tabThermo["enthalpy"].setQuantity(m_enthalpy_tab, "J/kmol");
     tabThermo["entropy"].setQuantity(m_entropy_tab, "J/kmol/K");
+    tabThermo["molar-volume"].setQuantity(m_molar_volume_tab, "m^3/kmol");
     phaseNode["tabulated-thermo"] = std::move(tabThermo);
 }
 
@@ -207,15 +199,18 @@ void BinarySolutionTabulatedThermo::diff(const vector_fp& inputData,
 
 void BinarySolutionTabulatedThermo::getPartialMolarVolumes(double* vbar) const
 {
-    vbar[m_kk_tab] = interpolate(moleFraction(m_kk_tab), m_partial_molar_volume_1_tab);
-    vbar[1-m_kk_tab] = interpolate(moleFraction(m_kk_tab),
-                                   m_partial_molar_volume_2_tab);
+    std::copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vbar);
 }
 
 void BinarySolutionTabulatedThermo::calcDensity()
 {
-    double vmol = interpolate(moleFraction(m_kk_tab), m_molar_volume_tab);
-    double dens = meanMolecularWeight()/vmol;
+    double Xtab = moleFraction(m_kk_tab);
+    double Vm = interpolate(Xtab, m_molar_volume_tab);
+    double dVdX_tab = interpolate(Xtab, m_derived_molar_volume_tab);
+    m_speciesMolarVolume[m_kk_tab] = Vm + (1 - Xtab) * dVdX_tab;
+    m_speciesMolarVolume[1-m_kk_tab] = Vm - Xtab * dVdX_tab;
+
+    double dens = meanMolecularWeight() / Vm;
 
     // Set the density in the parent State object directly, by calling the
     // Phase::assignDensity() function.

src/thermo/IdealSolidSolnPhase.cpp

@@ -43,7 +43,8 @@ doublereal IdealSolidSolnPhase::entropy_mole() const
 
 doublereal IdealSolidSolnPhase::gibbs_mole() const
 {
-    return RT() * (mean_X(gibbs_RT_ref()) + sum_xlogx());
+    double Pv = (pressure() - m_Pref)/molarDensity();
+    return RT() * (mean_X(gibbs_RT_ref()) + sum_xlogx()) + Pv;
 }
 
 doublereal IdealSolidSolnPhase::cp_mole() const

test/data/consistency-cases.yaml

@@ -77,9 +77,6 @@ ideal-condensed-1:
   setup:
     file: thermo-models.yaml
     phase: IdealSolidSolnPhase
-    known-failures:
-      g_eq_h_minus_Ts: "Inconsistent result when P != 1 atm. See GitHub Issue #1301"
-      g_eq_sum_gk_Xk: "Inconsistent result when P != 1 atm. See GitHub Issue #1301"
   states:
   - {T: 300, P: 0.1 atm, X: {sp1: 1.0}}
   - {T: 300, P: 1 atm, X: {sp1: 0.6, sp2: 0.4}}
@@ -90,9 +87,6 @@ ideal-condensed-2:
   setup:
     file: thermo-models.yaml
     phase: IdealSolidSolnPhase2
-    known-failures:
-      g_eq_h_minus_Ts: "Inconsistent result when P != 1 atm. See GitHub Issue #1301"
-      g_eq_sum_gk_Xk: "Inconsistent result when P != 1 atm. See GitHub Issue #1301"
   states:
   - {T: 300, P: 0.1 atm, X: {sp1: 1.0}}
   - {T: 400, P: 0.1 atm, X: {sp1: 0.01, sp2: 0.03, sp3: 0.94}}
@@ -103,12 +97,10 @@ binary-solution-tabulated:
     file: BinarySolutionTabulatedThermo.yaml
     phase: anode
     known-failures:
-      g_eq_h_minus_Ts: "Inconsistent results when P != 1 atm. See GitHub Issue #1301"
-      g_eq_sum_gk_Xk: "Inconsistent results when P != 1 atm. See GitHub Issue #1301"
-      v_eq_sum_vk_Xk: "Inconsistent results. See GitHub Issue #1302"
-      h_eq_sum_hk_Xk: "Inconsistent results. See GitHub Issue #1302"
-      s_eq_sum_sk_Xk: "Problems near composition limit. See GitHub Issue #1303"
-      gk_eq_hk_minus_T_sk: "Problems near composition limit. See GitHub Issue #1303"
+      g_eq_h_minus_Ts/4: "Problems near composition limit. See GitHub Issue #1303"
+      g_eq_sum_gk_Xk/4: "Problems near composition limit. See GitHub Issue #1303"
+      s_eq_sum_sk_Xk/4: "Problems near composition limit. See GitHub Issue #1303"
+      gk_eq_hk_minus_T_sk/[34]: "Problems near composition limit. See GitHub Issue #1303"
       gk0_eq_hk0_minus_T_sk0/[34]:
         "Problems near composition limit. See GitHub Issue #1303"
       chem_potentials_to_activities/3:

test/thermo/BinarySolutionTabulatedThermo_Test.cpp

@@ -58,15 +58,15 @@ TEST_F(BinarySolutionTabulatedThermo_Test,interp_s)
     test_phase->setState_TP(298.15, 101325.);
     // These expected results are purely a regression test
     const double expected_result[9] = {
-        3839.8896914480647,
-        5260.8983334513332,
-        5764.7097019695211,
-        7786.429533070881,
-        10411.474081913055,
-        15276.785945165157,
-        17900.243436157067,
-        22085.482962782506,
-        25989.144060372793
+        3839.8897013463629,
+        5260.8983501505654,
+        5764.7097249117705,
+        7786.4295622389736,
+        10411.474117668258,
+        15276.785988429341,
+        17900.243488318705,
+        22085.483025077028,
+        25989.144133134443
     };
 
     double xmin = 0.10;
@@ -90,15 +90,15 @@ TEST_F(BinarySolutionTabulatedThermo_Test,chem_potentials)
     test_phase->setState_TP(298.15,101325.);
     // These expected results are purely a regression test
     const double expected_result[9] = {
-        -19347891.714810669,
-        -14757822.388050893,
-        -12593133.605195494,
-        -12626837.865623865,
-        -12131010.479908356,
-        -10322881.86739888,
-        - 9573869.8636945337,
-        -10260863.826955771,
-        -10579827.307551134
+        -19347891.744322445,
+        -14757822.415520553,
+        -12593133.63125352,
+        -12626837.890922677,
+        -12131010.504991682,
+        -10322881.892878814,
+        -9573869.8901660107,
+        -10260863.854728648,
+        -10579827.336476315
     };
 
     double xmin = 0.10;
@@ -120,15 +120,15 @@ TEST_F(BinarySolutionTabulatedThermo_Test,partialMolarEntropies)
     test_phase->setState_TP(298.15,101325.);
     // These expected results are purely a regression test
     const double expected_result[9] = {
-        30514.752294683516,
-        21514.841983025333,
-        14848.02859501992,
-        15965.482659621264,
-        18272.567242414199,
-        24453.517437971925,
-        25299.003664716853,
-        28474.69918493319,
-        30810.094532734405
+        30514.752393666495,
+        21514.842075159024,
+        14848.028682418964,
+        15965.482744473897,
+        18272.567326544096,
+        24453.517523432041,
+        25299.003753502617,
+        28474.699278083881,
+        30810.094629749936
     };
 
     double xmin = 0.10;
@@ -176,15 +176,15 @@ TEST_F(BinarySolutionTabulatedThermo_Test,partialMolarVolumes)
     test_phase->setState_TP(298.15,101325.);
     // These expected results are purely a regression test
     const double expected_result[9] = {
-        0.041207972037360034,
-        0.038534004157808582,
-        0.036935982981359229,
-        0.036182506843878831,
-        0.035990796804076991,
-        0.036280986542177367,
-        0.036903215973399468,
-        0.037569211282710353,
-        0.038022737191326351
+        0.04120724363741,
+        0.03853288221791,
+        0.03693536558788,
+        0.03618236414389,
+        0.03599080437984,
+        0.03628136773515,
+        0.03690395850931,
+        0.03756972764230,
+        0.03802279519842
     };
 
     double xmin = 0.10;
@@ -196,7 +196,7 @@ TEST_F(BinarySolutionTabulatedThermo_Test,partialMolarVolumes)
     {
         set_defect_X(xmin + i*dx);
         test_phase->getPartialMolarVolumes(&partialMolarVolumes[0]);
-        EXPECT_NEAR(expected_result[i], partialMolarVolumes[0], 1.e-6);
+        EXPECT_NEAR(expected_result[i], partialMolarVolumes[0], 1.e-8);
     }
 }
 

test/thermo/phaseConstructors.cpp

@@ -446,7 +446,7 @@ TEST(IdealSolidSolnPhase, fromScratch)
     p.setState_TPX(500, 2e5, "sp1:0.1, sp2:0.89, sp3:0.01");
     EXPECT_NEAR(p.density(), 10.1787080, 1e-6);
     EXPECT_NEAR(p.enthalpy_mass(), -15642788.8547624, 1e-4);
-    EXPECT_NEAR(p.gibbs_mole(), -313642312.7114608, 1e-4);
+    EXPECT_NEAR(p.gibbs_mole(), -313513245.8114608, 1e-4);
 }
 
 static void set_hmw_interactions(HMWSoln& p) {

test/thermo/thermoFromYaml.cpp

@@ -390,7 +390,7 @@ TEST(ThermoFromYaml, IdealSolidSolnPhase)
     // Regression test following IdealSolidSolnPhase.fromScratch
     EXPECT_NEAR(thermo->density(), 10.1787080, 1e-6);
     EXPECT_NEAR(thermo->enthalpy_mass(), -15642788.8547624, 1e-4);
-    EXPECT_NEAR(thermo->gibbs_mole(), -313642312.7114608, 1e-4);
+    EXPECT_NEAR(thermo->gibbs_mole(), -313513245.8114608, 1e-4);
 
     // Test that molar enthalpy equals sum(h_k*X_k). Test first at default
     // pressure:

test/thermo/thermoToYaml.cpp

@@ -445,8 +445,8 @@ TEST_F(ThermoYamlRoundTrip, PengRobinson_crit_props)
 
 TEST_F(ThermoYamlRoundTrip, BinarySolutionTabulated)
 {
-    roundtrip("lithium_ion_battery.yaml", "cathode");
-    compareThermo(310, 2e5, "Li[cathode]:0.4, V[cathode]:0.6");
+    roundtrip("BinarySolutionTabulatedThermo.yaml", "anode");
+    compareThermo(310, 2e5, "Li[anode]:0.4, V[anode]:0.6");
 }
 
 TEST_F(ThermoYamlRoundTrip, Margules)