[Thermo] Add function 'resetHf298'

Using this instead of modifyHf298 to reset the thermo data to its original
state avoids round-off errors that otherwise make modifications to the
species thermo data irreversible.

include/cantera/thermo/ConstCpPoly.h

@@ -88,6 +88,7 @@ public:
 
     virtual doublereal reportHf298(doublereal* const h298 = 0) const;
     virtual void modifyOneHf298(const size_t k, const doublereal Hf298New);
+    virtual void resetHf298();
 
 protected:
     //! Base temperature
@@ -100,6 +101,8 @@ protected:
     doublereal m_s0_R;
     //! log of the t0 value
     doublereal m_logt0;
+    //! Original value of h0_R, restored by calling resetHf298()
+    double m_h0_R_orig;
 };
 
 }

include/cantera/thermo/GeneralSpeciesThermo.h

@@ -80,6 +80,7 @@ public:
     virtual doublereal reportOneHf298(const size_t k) const;
 
     virtual void modifyOneHf298(const size_t k, const doublereal Hf298New);
+    virtual void resetHf298(const size_t k);
 
 private:
     //! Provide the SpeciesthermoInterpType object

include/cantera/thermo/LatticeSolidPhase.h

@@ -436,6 +436,7 @@ public:
     void setLatticeMoleFractionsByName(int n, const std::string& x);
 
     virtual void modifyOneHf298SS(const size_t k, const doublereal Hf298New);
+    virtual void resetHf298(const size_t k=npos);
 
 protected:
     //! Current value of the pressure

include/cantera/thermo/NasaPoly1.h

@@ -59,6 +59,7 @@ public:
         : SpeciesThermoInterpType(tlow, thigh, pref)
         , m_coeff(coeffs, coeffs+7)
     {
+        m_coeff5_orig = m_coeff[5];
     }
 
     virtual SpeciesThermoInterpType*
@@ -163,9 +164,15 @@ public:
         m_coeff[5] += (delH) / GasConstant;
     }
 
+    virtual void resetHf298() {
+        m_coeff[5] = m_coeff5_orig;
+    }
+
 protected:
     //! array of polynomial coefficients, stored in the order [a0, ..., a6]
     vector_fp m_coeff;
+
+    double m_coeff5_orig;
 };
 
 }

include/cantera/thermo/NasaPoly2.h

@@ -136,6 +136,11 @@ public:
         return h;
     }
 
+    void resetHf298() {
+        mnp_low.resetHf298();
+        mnp_high.resetHf298();
+    }
+
     void modifyOneHf298(const size_t k, const doublereal Hf298New) {
         doublereal h298now = reportHf298(0);
         doublereal delH = Hf298New - h298now;

include/cantera/thermo/ShomatePoly.h

@@ -77,6 +77,7 @@ public:
         for (size_t i = 0; i < 7; i++) {
             m_coeff[i] = coeffs[i] * 1000 / GasConstant;
         }
+        m_coeff5_orig = m_coeff[5];
     }
 
     virtual SpeciesThermoInterpType*
@@ -168,9 +169,14 @@ public:
         m_coeff[5] += delH / (1e3 * GasConstant);
     }
 
+    virtual void resetHf298() {
+        m_coeff[5] = m_coeff5_orig;
+    }
+
 protected:
     //! Array of coefficients
     vector_fp m_coeff;
+    double m_coeff5_orig;
 };
 
 //! The Shomate polynomial parameterization for two temperature ranges for one
@@ -331,6 +337,11 @@ public:
         msp_high.modifyOneHf298(k, hnew);
     }
 
+    virtual void resetHf298() {
+        msp_low.resetHf298();
+        msp_high.resetHf298();
+    }
+
 protected:
     //! Midrange temperature (kelvin)
     doublereal m_midT;

include/cantera/thermo/SpeciesThermo.h

@@ -265,6 +265,13 @@ public:
      */
     virtual void modifyOneHf298(const size_t k, const doublereal Hf298New) = 0;
 
+    //! Restore the original heat of formation of one or more species
+    /*!
+     *  Resets changes made by modifyOneHf298(). If the species index is not
+     *  specified, the heats of formation for all species are restored.
+     */
+    virtual void resetHf298(const size_t k) = 0;
+
     //! Check if data for all species (0 through nSpecies-1) has been installed.
     bool ready(size_t nSpecies);
 

include/cantera/thermo/SpeciesThermoInterpType.h

@@ -13,6 +13,7 @@
 
 #include "cantera/base/ct_defs.h"
 #include "speciesThermoTypes.h"
+#include "cantera/base/ctexceptions.h"
 
 namespace Cantera
 {
@@ -278,6 +279,15 @@ public:
      */
     virtual void modifyOneHf298(const size_t k, const doublereal Hf298New);
 
+    //! Restore the original heat of formation for this species
+    /*!
+     *  Resets changes made by modifyOneHf298().
+     */
+    virtual void resetHf298() {
+        throw CanteraError("SpeciesThermoInterpType::resetHf298",
+                           "Not implemented");
+    }
+
 protected:
     //!  lowest valid temperature
     doublereal m_lowT;

include/cantera/thermo/ThermoPhase.h

@@ -178,6 +178,13 @@ public:
         invalidateCache();
     }
 
+    //! Restore the original heat of formation of one or more species
+    /*!
+     *  Resets changes made by modifyOneHf298SS(). If the species index is not
+     *  specified, the heats of formation for all species are restored.
+     */
+    virtual void resetHf298(const size_t k=npos);
+
     //! Maximum temperature for which the thermodynamic data for the species
     //! are valid.
     /*!

src/thermo/ConstCpPoly.cpp

@@ -28,6 +28,7 @@ ConstCpPoly::ConstCpPoly(double tlow, double thigh, double pref,
     m_s0_R = coeffs[2] / GasConstant;
     m_cp0_R = coeffs[3] / GasConstant;
     m_logt0 = log(m_t0);
+    m_h0_R_orig = m_h0_R;
 }
 
 SpeciesThermoInterpType*
@@ -103,4 +104,8 @@ void ConstCpPoly::modifyOneHf298(const size_t k, const doublereal Hf298New)
     m_h0_R += delH / GasConstant;
 }
 
+void ConstCpPoly::resetHf298() {
+    m_h0_R = m_h0_R_orig;
+}
+
 }

src/thermo/GeneralSpeciesThermo.cpp

@@ -254,4 +254,12 @@ void GeneralSpeciesThermo::modifyOneHf298(const size_t k, const doublereal Hf298
     }
 }
 
+void GeneralSpeciesThermo::resetHf298(const size_t k)
+{
+    SpeciesThermoInterpType* sp_ptr = provideSTIT(k);
+    if (sp_ptr) {
+        sp_ptr->resetHf298();
+    }
+}
+
 }

src/thermo/LatticeSolidPhase.cpp

@@ -459,4 +459,22 @@ void LatticeSolidPhase::modifyOneHf298SS(const size_t k, const doublereal Hf298N
     _updateThermo();
 }
 
+void LatticeSolidPhase::resetHf298(const size_t k)
+{
+    if (k != npos) {
+        for (size_t n = 0; n < m_lattice.size(); n++) {
+            if (lkstart_[n+1] < k) {
+                size_t kk = k-lkstart_[n];
+                m_lattice[n]->speciesThermo().resetHf298(kk);
+            }
+        }
+    } else {
+        for (size_t n = 0; n < m_lattice.size(); n++) {
+            m_lattice[n]->speciesThermo().resetHf298(npos);
+        }
+    }
+    invalidateCache();
+    _updateThermo();
+}
+
 } // End namespace Cantera

src/thermo/SpeciesThermoInterpType.cpp

@@ -6,7 +6,6 @@
 #include "cantera/thermo/SpeciesThermoInterpType.h"
 #include "cantera/thermo/VPSSMgr.h"
 #include "cantera/thermo/PDSS.h"
-#include "cantera/base/ctexceptions.h"
 
 namespace Cantera
 {

src/thermo/ThermoPhase.cpp

@@ -94,6 +94,17 @@ ThermoPhase* ThermoPhase::duplMyselfAsThermoPhase() const
     return new ThermoPhase(*this);
 }
 
+void ThermoPhase::resetHf298(size_t k) {
+    if (k != npos) {
+        m_spthermo->resetHf298(k);
+    } else {
+        for (size_t k = 0; k < nSpecies(); k++) {
+            m_spthermo->resetHf298(k);
+        }
+    }
+    invalidateCache();
+}
+
 int ThermoPhase::activityConvention() const
 {
     return cAC_CONVENTION_MOLAR;

test/thermo/thermoParameterizations.cpp

@@ -122,11 +122,15 @@ TEST(Shomate, modifyOneHf298)
 {
     ShomatePoly2 S(200, 6000, 101325, co2_shomate_coeffs);
 
-    EXPECT_NEAR(-393.5224e6, S.reportHf298(), 1e4);
+    double hf = S.reportHf298();
+    EXPECT_NEAR(-393.5224e6, hf, 1e4);
     double Htest = -400e6;
     S.modifyOneHf298(npos, Htest);
     double cp, h, s;
     S.updatePropertiesTemp(298.15, &cp, &h, &s);
     EXPECT_DOUBLE_EQ(Htest, h * 298.15 * GasConstant);
     EXPECT_DOUBLE_EQ(Htest, S.reportHf298());
+    S.resetHf298();
+    S.updatePropertiesTemp(298.15, &cp, &h, &s);
+    EXPECT_DOUBLE_EQ(hf, h * 298.15 * GasConstant);
 }