From e3afaf5e6151d85e1c8ff56a3a752fbc3b170955 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 2 Aug 2017 16:53:12 -0400 Subject: [PATCH] [Thermo] Make IonsFromNeutral constructible without XML --- .../cantera/thermo/IonsFromNeutralVPSSTP.h | 9 +- include/cantera/thermo/PDSS_IonsFromNeutral.h | 6 +- src/thermo/IonsFromNeutralVPSSTP.cpp | 311 +++++++----------- src/thermo/PDSS_IonsFromNeutral.cpp | 30 +- src/thermo/Species.cpp | 8 + test/thermo/phaseConstructors.cpp | 42 +++ 6 files changed, 193 insertions(+), 213 deletions(-) diff --git a/include/cantera/thermo/IonsFromNeutralVPSSTP.h b/include/cantera/thermo/IonsFromNeutralVPSSTP.h index 31522f9a0..bc5e3803b 100644 --- a/include/cantera/thermo/IonsFromNeutralVPSSTP.h +++ b/include/cantera/thermo/IonsFromNeutralVPSSTP.h @@ -271,18 +271,14 @@ public: //@} + virtual bool addSpecies(shared_ptr spec); void setNeutralMoleculePhase(shared_ptr neutral); shared_ptr getNeutralMoleculePhase(); virtual void initThermo(); - virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); virtual void setParametersFromXML(const XML_Node& thermoNode); private: - //! Initialize lengths of local variables after all species have - //! been identified. - void initLengths(); - //! Update the activity coefficients /*! * This function will be called to update the internally stored natural @@ -354,9 +350,6 @@ protected: //! Index of special species size_t indexSpecialSpecies_; - //! Index of special species - size_t indexSecondSpecialSpecies_; - //! Formula Matrix for composition of neutral molecules //! in terms of the molecules in this ThermoPhase /*! diff --git a/include/cantera/thermo/PDSS_IonsFromNeutral.h b/include/cantera/thermo/PDSS_IonsFromNeutral.h index 8181c80dd..e60b8636f 100644 --- a/include/cantera/thermo/PDSS_IonsFromNeutral.h +++ b/include/cantera/thermo/PDSS_IonsFromNeutral.h @@ -90,6 +90,8 @@ public: virtual bool useSTITbyPDSS() const { return true; } + void setNeutralSpeciesMultiplier(const std::string& species, double mult); + void setSpecialSpecies(bool special=true); void setParametersFromXML(const XML_Node& speciesNode); virtual void initThermo(); //@} @@ -100,7 +102,6 @@ protected: std::map neutralSpeciesMultipliers_; -public: //! Number of neutral molecule species that make up the stoichiometric //! vector for this species, in terms of calculating thermodynamic functions size_t numMult_; @@ -120,9 +121,6 @@ public: //! Vector of length equal to the number of species in the neutral molecule phase mutable vector_fp tmpNM; - - //! True if this species is the special species - int specialSpecies_; }; } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index 54fa9639f..3c0dfa7bb 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -30,7 +30,6 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() : ionSolnType_(cIonSolnType_SINGLEANION), numNeutralMoleculeSpecies_(0), indexSpecialSpecies_(npos), - indexSecondSpecialSpecies_(npos), neutralMoleculePhase_(0), geThermo(0) { @@ -40,8 +39,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile, const std::string& id_) : ionSolnType_(cIonSolnType_SINGLEANION), numNeutralMoleculeSpecies_(0), - indexSpecialSpecies_(npos), - indexSecondSpecialSpecies_(npos) + indexSpecialSpecies_(npos) { initThermoFile(inputFile, id_); } @@ -50,8 +48,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot, const std::string& id_) : ionSolnType_(cIonSolnType_SINGLEANION), numNeutralMoleculeSpecies_(0), - indexSpecialSpecies_(npos), - indexSecondSpecialSpecies_(npos) + indexSpecialSpecies_(npos) { importPhase(phaseRoot, this); } @@ -450,47 +447,6 @@ void IonsFromNeutralVPSSTP::compositionChanged() // ------------ Partial Molar Properties of the Solution ------------ -void IonsFromNeutralVPSSTP::initThermo() -{ - initLengths(); - GibbsExcessVPSSTP::initThermo(); -} - -void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr neutral) -{ - neutralMoleculePhase_ = neutral; - geThermo = dynamic_cast(neutralMoleculePhase_.get()); - numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies(); - fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk); - NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_); - muNeutralMolecule_.resize(numNeutralMoleculeSpecies_); - lnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); - dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); - dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); - dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); - dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0); - y_.resize(numNeutralMoleculeSpecies_, 0.0); - dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0); - dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0); -} - -shared_ptr IonsFromNeutralVPSSTP::getNeutralMoleculePhase() -{ - return neutralMoleculePhase_; -} - -void IonsFromNeutralVPSSTP::initLengths() -{ - moleFractions_.resize(m_kk); - fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk); - fm_invert_ionForNeutral.resize(m_kk); - cationList_.resize(m_kk); - anionList_.resize(m_kk); - passThroughList_.resize(m_kk); - moleFractionsTmp_.resize(m_kk); - m_work.resize(m_kk); -} - //! Return the factor overlap /*! * @param elnamesVN @@ -524,6 +480,129 @@ static double factorOverlap(const std::vector& elnamesVN , return fMax; } +void IonsFromNeutralVPSSTP::initThermo() +{ + size_t nElementsN = neutralMoleculePhase_->nElements(); + const std::vector& elnamesVN = neutralMoleculePhase_->elementNames(); + vector_fp elemVectorN(nElementsN); + + size_t nElementsI = nElements(); + const std::vector& elnamesVI = elementNames(); + vector_fp elemVectorI(nElementsI); + + for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + for (size_t m = 0; m < nElementsN; m++) { + elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m); + } + + for (size_t m = 0; m < nElementsI; m++) { + elemVectorI[m] = nAtoms(indexSpecialSpecies_, m); + } + double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, + elnamesVI ,elemVectorI, nElementsI); + if (fac > 0.0) { + for (size_t m = 0; m < nElementsN; m++) { + for (size_t mi = 0; mi < nElementsI; mi++) { + if (elnamesVN[m] == elnamesVI[mi]) { + elemVectorN[m] -= fac * elemVectorI[mi]; + } + + } + } + } + fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac; + + for (size_t k = 0; k < m_kk; k++) { + for (size_t m = 0; m < nElementsI; m++) { + elemVectorI[m] = nAtoms(k, m); + } + fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, + elnamesVI ,elemVectorI, nElementsI); + if (fac > 0.0) { + for (size_t m = 0; m < nElementsN; m++) { + for (size_t mi = 0; mi < nElementsI; mi++) { + if (elnamesVN[m] == elnamesVI[mi]) { + elemVectorN[m] -= fac * elemVectorI[mi]; + } + } + } + bool notTaken = true; + for (size_t iNeut = 0; iNeut < jNeut; iNeut++) { + if (fm_invert_ionForNeutral[k] == iNeut) { + notTaken = false; + } + } + if (notTaken) { + fm_invert_ionForNeutral[k] = jNeut; + } else { + throw CanteraError("IonsFromNeutralVPSSTP::initThermo", + "Simple formula matrix generation failed, one cation is shared between two salts"); + } + } + fm_neutralMolec_ions_[k + jNeut * m_kk] += fac; + } + + // Ok check the work + for (size_t m = 0; m < nElementsN; m++) { + if (fabs(elemVectorN[m]) > 1.0E-13) { + throw CanteraError("IonsFromNeutralVPSSTP::initThermo", + "Simple formula matrix generation failed"); + } + } + } + + GibbsExcessVPSSTP::initThermo(); +} + +void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr neutral) +{ + neutralMoleculePhase_ = neutral; + geThermo = dynamic_cast(neutralMoleculePhase_.get()); + numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies(); + fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk); + NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_); + muNeutralMolecule_.resize(numNeutralMoleculeSpecies_); + lnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0); + y_.resize(numNeutralMoleculeSpecies_, 0.0); + dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0); + dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0); +} + +shared_ptr IonsFromNeutralVPSSTP::getNeutralMoleculePhase() +{ + return neutralMoleculePhase_; +} + +bool IonsFromNeutralVPSSTP::addSpecies(shared_ptr spec) +{ + bool added = GibbsExcessVPSSTP::addSpecies(spec); + if (added) { + moleFractions_.push_back(0.0); + moleFractionsTmp_.push_back(0.0); + m_work.push_back(0.0); + fm_invert_ionForNeutral.push_back(npos); + fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk); + + if (spec->charge > 0) { + cationList_.push_back(m_kk-1); + } else if (spec->charge < 0) { + anionList_.push_back(m_kk-1); + } else { + passThroughList_.push_back(m_kk-1); + } + + if (spec->extra.hasKey("special_species") + && spec->extra["special_species"].asBool()) { + indexSpecialSpecies_ = m_kk - 1; + } + } + return added; +} + void IonsFromNeutralVPSSTP::setParametersFromXML(const XML_Node& thermoNode) { GibbsExcessVPSSTP::setParametersFromXML(thermoNode); @@ -543,146 +622,6 @@ void IonsFromNeutralVPSSTP::setParametersFromXML(const XML_Node& thermoNode) setNeutralMoleculePhase(shared_ptr(newPhase(*neut_ptr))); } -void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) -{ - if (id_.size() > 0 && phaseNode.id() != id_) { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "phasenode and Id are incompatible"); - } - - // Find the Thermo XML node - if (!phaseNode.hasChild("thermo")) { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "no thermo XML node"); - } - XML_Node& thermoNode = phaseNode.child("thermo"); - - // Make sure that the thermo model is IonsFromNeutralMolecule - if (!ba::iequals(thermoNode["model"], "ionsfromneutralmolecule")) { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "model name isn't IonsFromNeutralMolecule: " - + thermoNode["model"]); - } - - cationList_.clear(); - for (size_t k = 0; k < m_kk; k++) { - if (charge(k) > 0) { - cationList_.push_back(k); - } - } - - anionList_.clear(); - for (size_t k = 0; k < m_kk; k++) { - if (charge(k) < 0) { - anionList_.push_back(k); - } - } - - passThroughList_.clear(); - for (size_t k = 0; k < m_kk; k++) { - if (charge(k) == 0) { - passThroughList_.push_back(k); - } - } - - indexSpecialSpecies_ = npos; - for (size_t k = 0; k < m_kk; k++) { - PDSS_IonsFromNeutral* speciesSS = - dynamic_cast(providePDSS(k)); - if (!speciesSS) { - throw CanteraError("initThermoXML", "Dynamic cast failed"); - } - if (speciesSS->specialSpecies_ == 1) { - indexSpecialSpecies_ = k; - } - if (speciesSS->specialSpecies_ == 2) { - indexSecondSpecialSpecies_ = k; - } - } - - size_t nElementsN = neutralMoleculePhase_->nElements(); - const std::vector& elnamesVN = neutralMoleculePhase_->elementNames(); - vector_fp elemVectorN(nElementsN); - vector_fp elemVectorN_orig(nElementsN); - - size_t nElementsI = nElements(); - const std::vector& elnamesVI = elementNames(); - vector_fp elemVectorI(nElementsI); - - vector_fp fm_tmp(m_kk); - for (size_t k = 0; k < m_kk; k++) { - fm_invert_ionForNeutral[k] = npos; - } - for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { - for (size_t m = 0; m < nElementsN; m++) { - elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m); - } - elemVectorN_orig = elemVectorN; - fm_tmp.assign(m_kk, 0.0); - - for (size_t m = 0; m < nElementsI; m++) { - elemVectorI[m] = nAtoms(indexSpecialSpecies_, m); - } - double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, - elnamesVI ,elemVectorI, nElementsI); - if (fac > 0.0) { - for (size_t m = 0; m < nElementsN; m++) { - std::string mName = elnamesVN[m]; - for (size_t mi = 0; mi < nElementsI; mi++) { - std::string eName = elnamesVI[mi]; - if (mName == eName) { - elemVectorN[m] -= fac * elemVectorI[mi]; - } - - } - } - } - fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac; - - for (size_t k = 0; k < m_kk; k++) { - for (size_t m = 0; m < nElementsI; m++) { - elemVectorI[m] = nAtoms(k, m); - } - fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, - elnamesVI ,elemVectorI, nElementsI); - if (fac > 0.0) { - for (size_t m = 0; m < nElementsN; m++) { - std::string mName = elnamesVN[m]; - for (size_t mi = 0; mi < nElementsI; mi++) { - std::string eName = elnamesVI[mi]; - if (mName == eName) { - elemVectorN[m] -= fac * elemVectorI[mi]; - } - } - } - bool notTaken = true; - for (size_t iNeut = 0; iNeut < jNeut; iNeut++) { - if (fm_invert_ionForNeutral[k] == iNeut) { - notTaken = false; - } - } - if (notTaken) { - fm_invert_ionForNeutral[k] = jNeut; - } else { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "Simple formula matrix generation failed, one cation is shared between two salts"); - } - } - fm_neutralMolec_ions_[k + jNeut * m_kk] += fac; - } - - // Ok check the work - for (size_t m = 0; m < nElementsN; m++) { - if (fabs(elemVectorN[m]) > 1.0E-13) { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "Simple formula matrix generation failed"); - } - } - } - // This includes the setStateFromXML calls - GibbsExcessVPSSTP::initThermoXML(phaseNode, id_); -} - void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const { size_t icat, jNeut; diff --git a/src/thermo/PDSS_IonsFromNeutral.cpp b/src/thermo/PDSS_IonsFromNeutral.cpp index 489cc1c07..35226ae3b 100644 --- a/src/thermo/PDSS_IonsFromNeutral.cpp +++ b/src/thermo/PDSS_IonsFromNeutral.cpp @@ -21,7 +21,6 @@ PDSS_IonsFromNeutral::PDSS_IonsFromNeutral() : neutralMoleculePhase_(0) , numMult_(0) , add2RTln2_(true) - , specialSpecies_(0) { } @@ -30,6 +29,16 @@ void PDSS_IonsFromNeutral::setParent(VPStandardStateTP* phase, size_t k) neutralMoleculePhase_ = dynamic_cast(*phase).getNeutralMoleculePhase(); } +void PDSS_IonsFromNeutral::setNeutralSpeciesMultiplier(const std::string& species, double mult) +{ + neutralSpeciesMultipliers_[species] = mult; + numMult_++; +} + +void PDSS_IonsFromNeutral::setSpecialSpecies(bool special) { + add2RTln2_ = !special; +} + void PDSS_IonsFromNeutral::setParametersFromXML(const XML_Node& speciesNode) { PDSS::setParametersFromXML(speciesNode); @@ -49,21 +58,12 @@ void PDSS_IonsFromNeutral::setParametersFromXML(const XML_Node& speciesNode) "no Thermo::neutralSpeciesMultipliers Node for species " + speciesNode.name()); } - neutralSpeciesMultipliers_ = parseCompString(nsm->value()); - numMult_ = neutralSpeciesMultipliers_.size(); + for (auto& species_mult : parseCompString(nsm->value())) { + setNeutralSpeciesMultiplier(species_mult.first, species_mult.second); + } - specialSpecies_ = 0; - const XML_Node* ss = tn->findByName("specialSpecies"); - if (ss) { - specialSpecies_ = 1; - } - const XML_Node* sss = tn->findByName("secondSpecialSpecies"); - if (sss) { - specialSpecies_ = 2; - } - add2RTln2_ = true; - if (specialSpecies_ == 1) { - add2RTln2_ = false; + if (tn->findByName("specialSpecies")) { + setSpecialSpecies(); } } diff --git a/src/thermo/Species.cpp b/src/thermo/Species.cpp index ebfe3322d..0275cc4ea 100644 --- a/src/thermo/Species.cpp +++ b/src/thermo/Species.cpp @@ -66,6 +66,14 @@ shared_ptr newSpecies(const XML_Node& species_node) s->extra["molar_volume"] = getFloat(*stdstate, "molarVolume", "toSI"); } + // Extra data possibly used by IonsFromNeutralVPSSTP + const XML_Node* thermo = species_node.findByName("thermo"); + if (thermo && thermo->attrib("model") == "IonFromNeutral") { + if (thermo->hasChild("specialSpecies")) { + s->extra["special_species"] = true; + } + } + return s; } diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index 3cb539141..e06b8f6de 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -8,6 +8,7 @@ #include "cantera/thermo/WaterSSTP.h" #include "cantera/thermo/RedlichKwongMFTP.h" #include "cantera/thermo/IonsFromNeutralVPSSTP.h" +#include "cantera/thermo/PDSS_IonsFromNeutral.h" #include "cantera/thermo/IdealSolnGasVPSS.h" #include "cantera/thermo/IdealMolalSoln.h" #include "cantera/thermo/DebyeHuckel.h" @@ -112,6 +113,47 @@ TEST(IonsFromNeutralConstructor, fromXML) EXPECT_NEAR(mu[1], -2.88157298e+06, 1e-1); } +TEST(IonsFromNeutralConstructor, fromScratch) +{ + auto neutral = make_shared(); + neutral->addUndefinedElements(); + auto sKCl = make_shomate_species("KCl(L)", "K:1 Cl:1", kcl_shomate_coeffs); + neutral->addSpecies(sKCl); + std::unique_ptr ssKCl(new PDSS_ConstVol()); + ssKCl->setMolarVolume(0.03757); + neutral->installPDSS(0, std::move(ssKCl)); + neutral->initThermo(); + + IonsFromNeutralVPSSTP p; + p.addUndefinedElements(); + p.setNeutralMoleculePhase(neutral); + + auto sKp = make_shared("K+", parseCompString("K:1"), 1); + auto sClm = make_shared("Cl-", parseCompString("Cl:1"), -1); + sClm->extra["special_species"] = true; + p.addSpecies(sKp); + p.addSpecies(sClm); + std::unique_ptr ssKp(new PDSS_IonsFromNeutral()); + std::unique_ptr ssClm(new PDSS_IonsFromNeutral()); + ssKp->setNeutralSpeciesMultiplier("KCl(L)", 1.2); + ssClm->setNeutralSpeciesMultiplier("KCl(L)", 1.5); + ssClm->setSpecialSpecies(); + p.installPDSS(0, std::move(ssKp)); + p.installPDSS(1, std::move(ssClm)); + p.initThermo(); + + ASSERT_EQ((int) p.nSpecies(), 2); + p.setState_TPX(500, 2e5, "K+:0.1, Cl-:0.1"); + vector_fp mu(p.nSpecies()); + p.getChemPotentials(mu.data()); + + // Values for regression testing only -- same as XML test + EXPECT_NEAR(p.density(), 1984.3225978174073, 1e-6); + EXPECT_NEAR(p.enthalpy_mass(), -8035317241137.971, 1e-1); + EXPECT_NEAR(mu[0], -4.66404010e+08, 1e1); + EXPECT_NEAR(mu[1], -2.88157298e+06, 1e-1); +} + #ifndef HAS_NO_PYTHON // skip these tests if the Python converter is unavailable class CtiConversionTest : public testing::Test {