[Thermo] Allow instantiation of RedlichKisterVPSSTP without XML

include/cantera/thermo/RedlichKisterVPSSTP.h

@@ -357,6 +357,19 @@ public:
     virtual void initThermo();
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
+    //! Add a binary species interaction with the specified parameters
+    /*!
+     * @param speciesA         name of the first species
+     * @param speciesB         name of the second species
+     * @param excess_enthalpy  coefficients of the excess enthalpy polynomial
+     * @param n_enthalpy       number of excess enthalpy polynomial coefficients
+     * @param excess_entropy   coefficients of the excess entropy polynomial
+     * @param n_entropy        number of excess entropy polynomial coefficients
+     */
+    void addBinaryInteraction(const std::string& speciesA, const std::string& speciesB,
+        const double* excess_enthalpy, size_t n_enthalpy,
+        const double* excess_entropy, size_t n_entropy);
+
     //! @}
     //! @name  Derivatives of Thermodynamic Variables needed for Applications
     //! @{
@@ -381,13 +394,6 @@ private:
      */
     void readXMLBinarySpecies(XML_Node& xmlBinarySpecies);
 
-    //! Resize internal arrays within the object that depend upon the number
-    //! of binary Redlich-Kister interaction terms
-    /*!
-     *  @param num Number of binary Redlich-Kister interaction terms
-     */
-    void resizeNumInteractions(const size_t num);
-
     //! Initialize lengths of local variables after all species have been
     //! identified.
     void initLengths();

src/thermo/RedlichKisterVPSSTP.cpp

@@ -493,17 +493,6 @@ void RedlichKisterVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnAct
     }
 }
 
-void RedlichKisterVPSSTP::resizeNumInteractions(const size_t num)
-{
-    numBinaryInteractions_ = num;
-    m_pSpecies_A_ij.resize(num, npos);
-    m_pSpecies_B_ij.resize(num, npos);
-    m_N_ij.resize(num, npos);
-    m_HE_m_ij.resize(num);
-    m_SE_m_ij.resize(num);
-    dlnActCoeff_dX_.resize(num, num, 0.0);
-}
-
 void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
 {
     std::string xname = xmLBinarySpecies.name();
@@ -511,7 +500,6 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
         throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies",
                            "Incorrect name for processing this routine: " + xname);
     }
-    size_t Npoly = 0;
     vector_fp hParams, sParams;
     std::string iName = xmLBinarySpecies.attrib("speciesA");
     if (iName == "") {
@@ -529,50 +517,59 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
     if (iSpecies == npos) {
         return;
     }
-    string ispName = speciesName(iSpecies);
-    if (charge(iSpecies) != 0) {
-        throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "speciesA charge problem");
-    }
     size_t jSpecies = speciesIndex(jName);
     if (jSpecies == npos) {
         return;
     }
-    std::string jspName = speciesName(jSpecies);
-    if (charge(jSpecies) != 0) {
-        throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "speciesB charge problem");
-    }
 
     // Ok we have found a valid interaction
-    numBinaryInteractions_++;
-    size_t iSpot = numBinaryInteractions_ - 1;
-    m_pSpecies_A_ij.resize(numBinaryInteractions_);
-    m_pSpecies_B_ij.resize(numBinaryInteractions_);
-    m_pSpecies_A_ij[iSpot] = iSpecies;
-    m_pSpecies_B_ij[iSpot] = jSpecies;
-
     for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
         XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
         string nodeName = ba::to_lower_copy(xmlChild.name());
 
         // Process the binary species interaction child elements
         if (nodeName == "excessenthalpy") {
-            // Get the string containing all of the values
             getFloatArray(xmlChild, hParams, true, "toSI", "excessEnthalpy");
-            Npoly = std::max(hParams.size(), Npoly);
+        } else if (nodeName == "excessentropy") {
-        }
-
-        if (nodeName == "excessentropy") {
-            // Get the string containing all of the values
             getFloatArray(xmlChild, sParams, true, "toSI", "excessEntropy");
-            Npoly = std::max(sParams.size(), Npoly);
         }
     }
-    hParams.resize(Npoly, 0.0);
+    addBinaryInteraction(iName, jName, hParams.data(), hParams.size(),
-    sParams.resize(Npoly, 0.0);
+                         sParams.data(), sParams.size());
-    m_HE_m_ij.push_back(hParams);
+}
-    m_SE_m_ij.push_back(sParams);
+
-    m_N_ij.push_back(Npoly);
+void RedlichKisterVPSSTP::addBinaryInteraction(
-    resizeNumInteractions(numBinaryInteractions_);
+    const std::string& speciesA, const std::string& speciesB,
+    const double* excess_enthalpy, size_t n_enthalpy,
+    const double* excess_entropy, size_t n_entropy)
+{
+    size_t kA = speciesIndex(speciesA);
+    size_t kB = speciesIndex(speciesB);
+    if (kA == npos) {
+        throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
+            "Species '{}' not present in phase", speciesA);
+    } else if (kB == npos) {
+        throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
+            "Species '{}' not present in phase", speciesB);
+    }
+    if (charge(kA) != 0) {
+        throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
+            "Species '{}' should be neutral", speciesA);
+    } else if (charge(kB) != 0) {
+        throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
+            "Species '{}' should be neutral", speciesB);
+    }
+
+    m_pSpecies_A_ij.push_back(kA);
+    m_pSpecies_B_ij.push_back(kB);
+    m_HE_m_ij.emplace_back(excess_enthalpy, excess_enthalpy + n_enthalpy);
+    m_SE_m_ij.emplace_back(excess_entropy, excess_entropy + n_entropy);
+    size_t N = max(n_enthalpy, n_entropy);
+    m_HE_m_ij.back().resize(N, 0.0);
+    m_SE_m_ij.back().resize(N, 0.0);
+    m_N_ij.push_back(N);
+    dlnActCoeff_dX_.resize(N, N, 0.0);
+    numBinaryInteractions_++;
 }
 
 }

test/thermo/RedlichKisterTest.cpp

@@ -1,6 +1,9 @@
 #include "gtest/gtest.h"
 #include "cantera/thermo/RedlichKisterVPSSTP.h"
 #include "cantera/thermo/ThermoFactory.h"
+#include "cantera/thermo/ConstCpPoly.h"
+#include "cantera/base/stringUtils.h"
+#include "cantera/thermo/PDSS_IdealGas.h"
 
 namespace Cantera
 {
@@ -8,7 +11,9 @@ namespace Cantera
 class RedlichKister_Test : public testing::Test
 {
 public:
-    RedlichKister_Test() {
+    RedlichKister_Test() {}
+
+    void initXML() {
         test_phase.reset(newPhase("../data/RedlichKisterVPSSTP_valid.xml"));
     }
 
@@ -20,19 +25,8 @@ public:
     }
 
     std::unique_ptr<ThermoPhase> test_phase;
-};
 
-TEST_F(RedlichKister_Test, construct_from_xml)
+    const double expected_chempot[9] = {
-{
-    RedlichKisterVPSSTP* redlich_kister_phase = dynamic_cast<RedlichKisterVPSSTP*>(test_phase.get());
-    EXPECT_TRUE(redlich_kister_phase != NULL);
-}
-
-TEST_F(RedlichKister_Test, chem_potentials)
-{
-    test_phase->setState_TP(298.15, 101325.);
-
-    const double expected_result[9] = {
         -1.2791500420236044e+007,
         -1.2618554504124604e+007,
         -1.2445418272766629e+007,
@@ -44,6 +38,20 @@ TEST_F(RedlichKister_Test, chem_potentials)
         -1.1730895987035934e+007
     };
 
+};
+
+TEST_F(RedlichKister_Test, construct_from_xml)
+{
+    initXML();
+    RedlichKisterVPSSTP* redlich_kister_phase = dynamic_cast<RedlichKisterVPSSTP*>(test_phase.get());
+    ASSERT_TRUE(redlich_kister_phase != NULL);
+}
+
+TEST_F(RedlichKister_Test, chem_potentials)
+{
+    initXML();
+    test_phase->setState_TP(298.15, 101325.);
+
     double xmin = 0.6;
     double xmax = 0.9;
     int numSteps = 9;
@@ -53,12 +61,13 @@ TEST_F(RedlichKister_Test, chem_potentials)
     {
         set_r(xmin + i*dx);
         test_phase->getChemPotentials(&chemPotentials[0]);
-        EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6);
+        EXPECT_NEAR(expected_chempot[i], chemPotentials[0], 1.e-6);
     }
 }
 
 TEST_F(RedlichKister_Test, dlnActivities)
 {
+    initXML();
     test_phase->setState_TP(298.15, 101325.);
 
     const double expected_result[9] = {
@@ -89,6 +98,7 @@ TEST_F(RedlichKister_Test, dlnActivities)
 
 TEST_F(RedlichKister_Test, activityCoeffs)
 {
+    initXML();
     test_phase->setState_TP(298., 1.);
 
     // Test that mu0 + RT log(activityCoeff * MoleFrac) == mu
@@ -115,12 +125,14 @@ TEST_F(RedlichKister_Test, activityCoeffs)
 
 TEST_F(RedlichKister_Test, standardConcentrations)
 {
+    initXML();
     EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(0));
     EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(1));
 }
 
 TEST_F(RedlichKister_Test, activityConcentrations)
 {
+    initXML();
     // Check to make sure activityConcentration_i == standardConcentration_i * gamma_i * X_i
     vector_fp standardConcs(2);
     vector_fp activityCoeffs(2);
@@ -144,4 +156,48 @@ TEST_F(RedlichKister_Test, activityConcentrations)
     }
 }
 
+TEST_F(RedlichKister_Test, fromScratch)
+{
+    test_phase.reset(new RedlichKisterVPSSTP());
+    RedlichKisterVPSSTP& rk = dynamic_cast<RedlichKisterVPSSTP&>(*test_phase);
+
+    auto sLiC6 = make_shared<Species>("Li(C6)", parseCompString("C:6 Li:1"));
+    double coeffs1[] = {298.15, -11.65e6, 0.0, 0.0};
+    sLiC6->thermo.reset(new ConstCpPoly(100, 5000, 101325, coeffs1));
+
+    auto sVC6 = make_shared<Species>("V(C6)", parseCompString("C:6"));
+    double coeffs2[] = {298.15, 0.0, 0.0, 0.0};
+    sVC6->thermo.reset(new ConstCpPoly(250, 800, 101325, coeffs2));
+
+    rk.addUndefinedElements();
+    rk.addSpecies(sLiC6);
+    rk.addSpecies(sVC6);
+
+    std::unique_ptr<PDSS> ssLiC6(new PDSS_IdealGas());
+    rk.installPDSS(0, std::move(ssLiC6));
+
+    std::unique_ptr<PDSS> ssVC6(new PDSS_IdealGas());
+    rk.installPDSS(1, std::move(ssVC6));
+
+    double hcoeffs[] = {-3.268E6, 3.955E6, -4.573E6, 6.147E6, -3.339E6, 1.117E7,
+                        2.997E5, -4.866E7, 1.362E5, 1.373E8, -2.129E7, -1.722E8,
+                        3.956E7, 9.302E7, -3.280E7};
+    double scoeffs[] = {0.0};
+
+    rk.addBinaryInteraction("Li(C6)", "V(C6)", hcoeffs, 15, scoeffs, 1);
+    rk.initThermo();
+    rk.setState_TPX(298.15, 101325, "Li(C6):0.6,V(C6):0.4");
+
+    double xmin = 0.6;
+    double xmax = 0.9;
+    int numSteps = 9;
+    vector_fp chemPotentials(2);
+    for(int i=0; i < 9; ++i)
+    {
+        set_r(xmin + i*(xmax-xmin)/(numSteps-1));
+        test_phase->getChemPotentials(&chemPotentials[0]);
+        EXPECT_NEAR(expected_chempot[i], chemPotentials[0], 1.e-6);
+    }
+}
+
 };