From f8ef4a8b2b3efc7392b5cdfd2a378adddf7579b8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 19 Feb 2017 19:16:27 -0500 Subject: [PATCH] [Thermo] Allow instantiation of RedlichKisterVPSSTP without XML --- include/cantera/thermo/RedlichKisterVPSSTP.h | 20 +++-- src/thermo/RedlichKisterVPSSTP.cpp | 77 +++++++++--------- test/thermo/RedlichKisterTest.cpp | 84 ++++++++++++++++---- 3 files changed, 120 insertions(+), 61 deletions(-) diff --git a/include/cantera/thermo/RedlichKisterVPSSTP.h b/include/cantera/thermo/RedlichKisterVPSSTP.h index b518a1a53..4f77d53bc 100644 --- a/include/cantera/thermo/RedlichKisterVPSSTP.h +++ b/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(); diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp index 8d352055d..1354e0ae3 100644 --- a/src/thermo/RedlichKisterVPSSTP.cpp +++ b/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); - } - - if (nodeName == "excessentropy") { - // Get the string containing all of the values + } else if (nodeName == "excessentropy") { getFloatArray(xmlChild, sParams, true, "toSI", "excessEntropy"); - Npoly = std::max(sParams.size(), Npoly); } } - hParams.resize(Npoly, 0.0); - sParams.resize(Npoly, 0.0); - m_HE_m_ij.push_back(hParams); - m_SE_m_ij.push_back(sParams); - m_N_ij.push_back(Npoly); - resizeNumInteractions(numBinaryInteractions_); + addBinaryInteraction(iName, jName, hParams.data(), hParams.size(), + sParams.data(), sParams.size()); +} + +void RedlichKisterVPSSTP::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) +{ + 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_++; } } diff --git a/test/thermo/RedlichKisterTest.cpp b/test/thermo/RedlichKisterTest.cpp index 6192cd01c..4543564f3 100644 --- a/test/thermo/RedlichKisterTest.cpp +++ b/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 test_phase; -}; -TEST_F(RedlichKister_Test, construct_from_xml) -{ - RedlichKisterVPSSTP* redlich_kister_phase = dynamic_cast(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] = { + const double expected_chempot[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(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(*test_phase); + + auto sLiC6 = make_shared("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("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 ssLiC6(new PDSS_IdealGas()); + rk.installPDSS(0, std::move(ssLiC6)); + + std::unique_ptr 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); + } +} + };