From 9c4a0baa558510bf0fe371a1e5eac27d3d5e2d5b Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 15 Apr 2016 20:18:09 -0400 Subject: [PATCH] [Thermo] Simplify adding species for most phase types Where possible, extend arrays as species are added rather than requiring a later call to initThermo(). For phases that do not require any data except that which is included in the Species objects themselves (notably, this includes IdealGasPhase), species can now be added dynamically without affecting the phase state. --- include/cantera/thermo/ConstDensityThermo.h | 2 +- include/cantera/thermo/DebyeHuckel.h | 17 +----- include/cantera/thermo/GibbsExcessVPSSTP.h | 7 +-- include/cantera/thermo/IdealGasPhase.h | 2 +- include/cantera/thermo/IdealMolalSoln.h | 8 +-- include/cantera/thermo/IdealSolidSolnPhase.h | 5 +- include/cantera/thermo/IdealSolnGasVPSS.h | 9 +--- include/cantera/thermo/LatticePhase.h | 2 +- include/cantera/thermo/LatticeSolidPhase.h | 4 +- include/cantera/thermo/MixtureFugacityTP.h | 8 +-- include/cantera/thermo/MolalityVPSSTP.h | 6 +-- include/cantera/thermo/Phase.h | 7 +++ include/cantera/thermo/RedlichKwongMFTP.h | 7 +-- include/cantera/thermo/SingleSpeciesTP.h | 8 +-- include/cantera/thermo/SurfPhase.h | 3 +- include/cantera/thermo/ThermoPhase.h | 1 - src/thermo/ConstDensityThermo.cpp | 17 +++--- src/thermo/DebyeHuckel.cpp | 55 ++++++++------------ src/thermo/GibbsExcessVPSSTP.cpp | 33 ++++++------ src/thermo/IdealGasPhase.cpp | 23 ++++---- src/thermo/IdealMolalSoln.cpp | 20 ++++--- src/thermo/IdealSolidSolnPhase.cpp | 45 ++++++++-------- src/thermo/IdealSolnGasVPSS.cpp | 16 +++--- src/thermo/LatticePhase.cpp | 22 ++++---- src/thermo/LatticeSolidPhase.cpp | 14 ++--- src/thermo/MixtureFugacityTP.cpp | 33 +++++------- src/thermo/MolalityVPSSTP.cpp | 27 ++++------ src/thermo/Phase.cpp | 1 + src/thermo/RedlichKwongMFTP.cpp | 41 +++++++-------- src/thermo/SingleSpeciesTP.cpp | 32 +++++------- src/thermo/SurfPhase.cpp | 27 +++++----- src/thermo/ThermoPhase.cpp | 7 +-- 32 files changed, 216 insertions(+), 293 deletions(-) diff --git a/include/cantera/thermo/ConstDensityThermo.h b/include/cantera/thermo/ConstDensityThermo.h index d89353bcb..5153e572a 100644 --- a/include/cantera/thermo/ConstDensityThermo.h +++ b/include/cantera/thermo/ConstDensityThermo.h @@ -121,7 +121,7 @@ public: return m_cp0_R; } - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); virtual void setParameters(int n, doublereal* const c) { setDensity(c[0]); diff --git a/include/cantera/thermo/DebyeHuckel.h b/include/cantera/thermo/DebyeHuckel.h index 869e3fca2..18e21aff6 100644 --- a/include/cantera/thermo/DebyeHuckel.h +++ b/include/cantera/thermo/DebyeHuckel.h @@ -829,7 +829,7 @@ public: * -------------- Utilities ------------------------------- */ - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); //! Return the Debye Huckel constant as a function of temperature @@ -1112,14 +1112,6 @@ protected: */ vector_fp m_B_Dot; - /** - * These are coefficients to describe the increase in activity coeff for - * non-polar molecules due to the electrolyte becoming stronger (the - * so-called salt-out effect) - */ - vector_fp m_npActCoeff; - - //! Pointer to the Water standard state object /*! * derived from the equation of state for water. @@ -1180,13 +1172,6 @@ protected: mutable vector_fp m_dlnActCoeffMolaldP; private: - //! Initialize the internal lengths. - /*! - * This internal function adjusts the lengths of arrays based on - * the number of species. - */ - void initLengths(); - //! Calculate the log activity coefficients /*! * This function updates the internally stored natural logarithm of the diff --git a/include/cantera/thermo/GibbsExcessVPSSTP.h b/include/cantera/thermo/GibbsExcessVPSSTP.h index 4bf61d420..4520275ac 100644 --- a/include/cantera/thermo/GibbsExcessVPSSTP.h +++ b/include/cantera/thermo/GibbsExcessVPSSTP.h @@ -244,12 +244,7 @@ public: virtual void setConcentrations(const doublereal* const c); //@} - virtual void initThermo(); - -private: - //! Initialize lengths of local variables after all species have - //! been identified. - void initLengths(); + virtual bool addSpecies(shared_ptr spec); protected: //! utility routine to check mole fraction sum diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index 23ac2e799..f3f8c1e04 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -605,7 +605,7 @@ public: //@} - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); virtual void setToEquilState(const doublereal* lambda_RT); protected: diff --git a/include/cantera/thermo/IdealMolalSoln.h b/include/cantera/thermo/IdealMolalSoln.h index 3e6906da9..b5f81ebd7 100644 --- a/include/cantera/thermo/IdealMolalSoln.h +++ b/include/cantera/thermo/IdealMolalSoln.h @@ -412,7 +412,7 @@ public: // -------------- Utilities ------------------------------- - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id=""); @@ -513,12 +513,6 @@ private: */ void s_updateIMS_lnMolalityActCoeff() const; - //! This internal function adjusts the lengths of arrays. - /*! - * This function is not virtual nor is it inherited - */ - void initLengths(); - //! Calculate parameters for cutoff treatments of activity coefficients /*! * Some cutoff treatments for the activity coefficients actually require diff --git a/include/cantera/thermo/IdealSolidSolnPhase.h b/include/cantera/thermo/IdealSolidSolnPhase.h index 8ef4eb28f..eaf269784 100644 --- a/include/cantera/thermo/IdealSolidSolnPhase.h +++ b/include/cantera/thermo/IdealSolidSolnPhase.h @@ -619,8 +619,8 @@ public: /// @name Utility Functions //@{ + virtual bool addSpecies(shared_ptr spec); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); - virtual void setToEquilState(const doublereal* lambda_RT); /** @@ -718,9 +718,6 @@ private: */ void _updateThermo() const; - //! This internal function adjusts the lengths of arrays - void initLengths(); - //@} }; } diff --git a/include/cantera/thermo/IdealSolnGasVPSS.h b/include/cantera/thermo/IdealSolnGasVPSS.h index bf0203d0d..c097e6525 100644 --- a/include/cantera/thermo/IdealSolnGasVPSS.h +++ b/include/cantera/thermo/IdealSolnGasVPSS.h @@ -144,18 +144,11 @@ public: */ //@{ + virtual bool addSpecies(shared_ptr spec); virtual void setParametersFromXML(const XML_Node& thermoNode); - virtual void initThermo(); virtual void setToEquilState(const doublereal* lambda_RT); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); -private: - //! @internal Initialize the internal lengths in this object. - /*! - * Note this is not a virtual function and only handles this object - */ - void initLengths(); - //@} protected: diff --git a/include/cantera/thermo/LatticePhase.h b/include/cantera/thermo/LatticePhase.h index bcd4f2928..e7f165cad 100644 --- a/include/cantera/thermo/LatticePhase.h +++ b/include/cantera/thermo/LatticePhase.h @@ -625,7 +625,7 @@ public: /// @name Utilities for Initialization of the Object //@{ - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); //! Set the equation of state parameters from the argument list diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h index f0bf0d56e..af237589b 100644 --- a/include/cantera/thermo/LatticeSolidPhase.h +++ b/include/cantera/thermo/LatticeSolidPhase.h @@ -421,11 +421,9 @@ public: virtual void getGibbs_RT_ref(doublereal* grt) const; virtual void getGibbs_ref(doublereal* g) const; + virtual bool addSpecies(shared_ptr spec); virtual void initThermo(); - //! Initialize vectors that depend on the number of species and sublattices - void initLengths(); - //! Add in species from Slave phases /*! * This hook is used for cSS_CONVENTION_SLAVE phases diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h index 2b4932f79..bdea15983 100644 --- a/include/cantera/thermo/MixtureFugacityTP.h +++ b/include/cantera/thermo/MixtureFugacityTP.h @@ -372,14 +372,8 @@ public: */ //@{ + virtual bool addSpecies(shared_ptr spec); virtual void setStateFromXML(const XML_Node& state); - virtual void initThermo(); - virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); - -private: - //! @internal Initialize the internal lengths in this object. - void initLengths(); - //@} protected: //! @name Special Functions for fugacity classes diff --git a/include/cantera/thermo/MolalityVPSSTP.h b/include/cantera/thermo/MolalityVPSSTP.h index 8a914fb27..6ad2d99c0 100644 --- a/include/cantera/thermo/MolalityVPSSTP.h +++ b/include/cantera/thermo/MolalityVPSSTP.h @@ -485,8 +485,8 @@ public: /// see importPhase(). //@{ + virtual bool addSpecies(shared_ptr spec); virtual void initThermo(); - virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); //@} @@ -567,10 +567,6 @@ private: */ virtual size_t findCLMIndex() const; - //! Initialize lengths of local variables after all species have been - //! identified. - void initLengths(); - protected: //! Index of the solvent. Currently the index of the solvent is hard-coded diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 96eb1827d..180cb8d13 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -686,6 +686,13 @@ public: //! Add a Species to this Phase. Returns `true` if the species was //! successfully added, or `false` if the species was ignored. + //! + //! Derived classes which need to size arrays according to the number of + //! species should overload this method. The derived class implementation + //! should call the base class method, and, if this returns `true` + //! (indicating that the species has been added), adjust their array sizes + //! accordingly. + //! //! @see ignoreUndefinedElements addUndefinedElements throwUndefinedElements virtual bool addSpecies(shared_ptr spec); diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h index 688660452..b7d3d09e3 100644 --- a/include/cantera/thermo/RedlichKwongMFTP.h +++ b/include/cantera/thermo/RedlichKwongMFTP.h @@ -190,8 +190,8 @@ public: */ //@{ + virtual bool addSpecies(shared_ptr spec); virtual void setParametersFromXML(const XML_Node& thermoNode); - virtual void initThermo(); virtual void setToEquilState(const doublereal* lambda_RT); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); @@ -211,11 +211,6 @@ private: */ void readXMLCrossFluid(XML_Node& pureFluidParam); - //! @internal Initialize the internal lengths in this object. - /*! - * Note this is not a virtual function and only handles this object - */ - void initLengths(); // @} protected: diff --git a/include/cantera/thermo/SingleSpeciesTP.h b/include/cantera/thermo/SingleSpeciesTP.h index 32ed6dd0f..c5519461e 100644 --- a/include/cantera/thermo/SingleSpeciesTP.h +++ b/include/cantera/thermo/SingleSpeciesTP.h @@ -252,13 +252,7 @@ public: doublereal tol = 1.e-8); //@} - /** - * @copydoc ThermoPhase::initThermo - * - * This version sets the mole fraction vector to x[0] = 1.0, and then - * calls the ThermoPhase::initThermo() function. - */ - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); protected: //! The current pressure of the solution (Pa). It gets initialized to 1 atm. diff --git a/include/cantera/thermo/SurfPhase.h b/include/cantera/thermo/SurfPhase.h index b5053b2dc..8c1c14ba0 100644 --- a/include/cantera/thermo/SurfPhase.h +++ b/include/cantera/thermo/SurfPhase.h @@ -293,8 +293,7 @@ public: */ virtual void setParametersFromXML(const XML_Node& thermoData); - //! Initialize the SurfPhase object after all species have been set up - virtual void initThermo(); + virtual bool addSpecies(shared_ptr spec); //! Set the initial state of the Surface Phase from an XML_Node /*! diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index ab726e3c3..e3f4eec9c 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -1381,7 +1381,6 @@ public: */ //@{ - using Phase::addSpecies; virtual bool addSpecies(shared_ptr spec); virtual void modifySpecies(size_t k, shared_ptr spec); diff --git a/src/thermo/ConstDensityThermo.cpp b/src/thermo/ConstDensityThermo.cpp index 96e9abff2..96d8a2aec 100644 --- a/src/thermo/ConstDensityThermo.cpp +++ b/src/thermo/ConstDensityThermo.cpp @@ -109,14 +109,17 @@ void ConstDensityThermo::getStandardChemPotentials(doublereal* mu0) const getPureGibbs(mu0); } -void ConstDensityThermo::initThermo() +bool ConstDensityThermo::addSpecies(shared_ptr spec) { - ThermoPhase::initThermo(); - m_h0_RT.resize(m_kk); - m_g0_RT.resize(m_kk); - m_cp0_R.resize(m_kk); - m_s0_R.resize(m_kk); - m_pp.resize(m_kk); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + m_h0_RT.push_back(0.0); + m_g0_RT.push_back(0.0); + m_cp0_R.push_back(0.0); + m_s0_R.push_back(0.0); + m_pp.push_back(0.0); + } + return added; } void ConstDensityThermo::_updateThermo() const diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index 7c283b52c..2e9eada3b 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -100,7 +100,6 @@ DebyeHuckel& DebyeHuckel::operator=(const DebyeHuckel& b) m_formDH = b.m_formDH; m_formGC = b.m_formGC; m_Aionic = b.m_Aionic; - m_npActCoeff = b.m_npActCoeff; m_IionicMolality = b.m_IionicMolality; m_maxIionicStrength = b.m_maxIionicStrength; m_useHelgesonFixedForm= b.m_useHelgesonFixedForm; @@ -109,7 +108,6 @@ DebyeHuckel& DebyeHuckel::operator=(const DebyeHuckel& b) m_A_Debye = b.m_A_Debye; m_B_Debye = b.m_B_Debye; m_B_Dot = b.m_B_Dot; - m_npActCoeff = b.m_npActCoeff; // This is an internal shallow copy of the PDSS_Water pointer m_waterSS = dynamic_cast(providePDSS(0)); @@ -400,16 +398,6 @@ void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const // -------------- Utilities ------------------------------- -void DebyeHuckel::initThermo() -{ - MolalityVPSSTP::initThermo(); - m_npActCoeff.resize(3); - m_npActCoeff[0] = 0.1127; - m_npActCoeff[1] = -0.01049; - m_npActCoeff[2] = 1.545E-3; - initLengths(); -} - //! Utility function to assign an integer value from a string for the //! ElectrolyteSpeciesType field. /*! @@ -542,10 +530,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_) " should be first species"); } - // Initialize all of the lengths of arrays in the object now that we know - // what species are in the phase. - initThermo(); - // Now go get the specification of the standard states for species in the // solution. This includes the molar volumes data blocks for incompressible // species. @@ -727,6 +711,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_) if (acNode.hasChild("DHBetaMatrix")) { if (m_formDH == DHFORM_BETAIJ || m_formDH == DHFORM_PITZER_BETAIJ) { + m_Beta_ij.resize(m_kk, m_kk, 0.0); XML_Node& irNode = acNode.child("DHBetaMatrix"); const vector& sn = speciesNames(); getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true); @@ -942,31 +927,35 @@ double DebyeHuckel::AionicRadius(int k) const // ------------ Private and Restricted Functions ------------------ -void DebyeHuckel::initLengths() +bool DebyeHuckel::addSpecies(shared_ptr spec) { - m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral); - m_speciesSize.resize(m_kk); - m_Aionic.resize(m_kk, 0.0); - m_lnActCoeffMolal.resize(m_kk, 0.0); - m_dlnActCoeffMolaldT.resize(m_kk, 0.0); - m_d2lnActCoeffMolaldT2.resize(m_kk, 0.0); - m_dlnActCoeffMolaldP.resize(m_kk, 0.0); - m_B_Dot.resize(m_kk, 0.0); - m_pp.resize(m_kk, 0.0); - m_tmpV.resize(m_kk, 0.0); - if (m_formDH == DHFORM_BETAIJ || - m_formDH == DHFORM_PITZER_BETAIJ) { - m_Beta_ij.resize(m_kk, m_kk, 0.0); + bool added = MolalityVPSSTP::addSpecies(spec); + if (added) { + m_electrolyteSpeciesType.push_back(cEST_polarNeutral); + m_speciesSize.push_back(0.0); + m_Aionic.push_back(0.0); + m_lnActCoeffMolal.push_back(0.0); + m_dlnActCoeffMolaldT.push_back(0.0); + m_d2lnActCoeffMolaldT2.push_back(0.0); + m_dlnActCoeffMolaldP.push_back(0.0); + m_B_Dot.push_back(0.0); + m_pp.push_back(0.0); + m_tmpV.push_back(0.0); } + return added; } double DebyeHuckel::_nonpolarActCoeff(double IionicMolality) const { + // These are coefficients to describe the increase in activity coeff for + // non-polar molecules due to the electrolyte becoming stronger (the so- + // called salt-out effect) + const static double npActCoeff[] = {0.1127, -0.01049, 1.545E-3}; double I2 = IionicMolality * IionicMolality; double l10actCoeff = - m_npActCoeff[0] * IionicMolality + - m_npActCoeff[1] * I2 + - m_npActCoeff[2] * I2 * IionicMolality; + npActCoeff[0] * IionicMolality + + npActCoeff[1] * I2 + + npActCoeff[2] * I2 * IionicMolality; return pow(10.0 , l10actCoeff); } diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index 1543a612d..683a17c03 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -157,23 +157,24 @@ double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const return norm; } -void GibbsExcessVPSSTP::initThermo() +bool GibbsExcessVPSSTP::addSpecies(shared_ptr spec) { - initLengths(); - VPStandardStateTP::initThermo(); - getMoleFractions(moleFractions_.data()); -} - -void GibbsExcessVPSSTP::initLengths() -{ - moleFractions_.resize(m_kk); - lnActCoeff_Scaled_.resize(m_kk); - dlnActCoeffdT_Scaled_.resize(m_kk); - d2lnActCoeffdT2_Scaled_.resize(m_kk); - dlnActCoeffdlnX_diag_.resize(m_kk); - dlnActCoeffdlnN_diag_.resize(m_kk); - dlnActCoeffdlnN_.resize(m_kk, m_kk); - m_pp.resize(m_kk); + bool added = VPStandardStateTP::addSpecies(spec); + if (added) { + if (m_kk == 1) { + moleFractions_.push_back(1.0); + } else { + moleFractions_.push_back(0.0); + } + lnActCoeff_Scaled_.push_back(0.0); + dlnActCoeffdT_Scaled_.push_back(0.0); + d2lnActCoeffdT2_Scaled_.push_back(0.0); + dlnActCoeffdlnX_diag_.push_back(0.0); + dlnActCoeffdlnN_diag_.push_back(0.0); + dlnActCoeffdlnN_.resize(m_kk, m_kk); + m_pp.push_back(0.0); + } + return added; } } // end of namespace Cantera diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index ce0ac3990..96a0a48ca 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -261,16 +261,21 @@ void IdealGasPhase::getStandardVolumes_ref(doublereal* vol) const } } -void IdealGasPhase::initThermo() +bool IdealGasPhase::addSpecies(shared_ptr spec) { - ThermoPhase::initThermo(); - m_p0 = refPressure(); - m_h0_RT.resize(m_kk); - m_g0_RT.resize(m_kk); - m_expg0_RT.resize(m_kk); - m_cp0_R.resize(m_kk); - m_s0_R.resize(m_kk); - m_pp.resize(m_kk); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + if (m_kk == 1) { + m_p0 = refPressure(); + } + m_h0_RT.push_back(0.0); + m_g0_RT.push_back(0.0); + m_expg0_RT.push_back(0.0); + m_cp0_R.push_back(0.0); + m_s0_R.push_back(0.0); + m_pp.push_back(0.0); + } + return added; } void IdealGasPhase::setToEquilState(const doublereal* mu_RT) diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp index 2a054f187..60e392c77 100644 --- a/src/thermo/IdealMolalSoln.cpp +++ b/src/thermo/IdealMolalSoln.cpp @@ -393,10 +393,16 @@ void IdealMolalSoln::getPartialMolarCp(doublereal* cpbar) const // -------------- Utilities ------------------------------- -void IdealMolalSoln::initThermo() +bool IdealMolalSoln::addSpecies(shared_ptr spec) { - initLengths(); - MolalityVPSSTP::initThermo(); + bool added = MolalityVPSSTP::addSpecies(spec); + if (added) { + m_pp.push_back(0.0); + m_speciesMolarVolume.push_back(0.0); + m_tmpV.push_back(0.0); + IMS_lnActCoeffMolal_.push_back(0.0); + } + return added; } void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) @@ -641,14 +647,6 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const } } -void IdealMolalSoln::initLengths() -{ - m_pp.resize(m_kk); - m_speciesMolarVolume.resize(m_kk); - m_tmpV.resize(m_kk); - IMS_lnActCoeffMolal_.resize(m_kk); -} - void IdealMolalSoln::calcIMSCutoffParams_() { IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_); diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp index 685a05718..060dfb485 100644 --- a/src/thermo/IdealSolidSolnPhase.cpp +++ b/src/thermo/IdealSolidSolnPhase.cpp @@ -450,6 +450,30 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const // Utility Functions +bool IdealSolidSolnPhase::addSpecies(shared_ptr spec) +{ + bool added = ThermoPhase::addSpecies(spec); + if (added) { + if (m_kk == 1) { + // Obtain the reference pressure by calling the ThermoPhase function + // refPressure, which in turn calls the species thermo reference + // pressure function of the same name. + m_Pref = refPressure(); + } + + m_h0_RT.push_back(0.0); + m_g0_RT.push_back(0.0); + m_expg0_RT.push_back(0.0); + m_cp0_R.push_back(0.0); + m_s0_R.push_back(0.0); + m_pe.push_back(0.0);; + m_pp.push_back(0.0); + m_speciesMolarVolume.push_back(0.0); + } + return added; +} + + void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_) { if (id_.size() > 0 && phaseNode.id() != id_) { @@ -495,10 +519,6 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& "Unspecified standardConc model"); } - // Initialize all of the lengths now that we know how many species - // there are in the phase. - initLengths(); - // Now go get the molar volumes XML_Node& speciesList = phaseNode.child("speciesArray"); XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], @@ -514,23 +534,6 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& ThermoPhase::initThermoXML(phaseNode, id_); } -void IdealSolidSolnPhase::initLengths() -{ - // Obtain the reference pressure by calling the ThermoPhase function - // refPressure, which in turn calls the species thermo reference pressure - // function of the same name. - m_Pref = refPressure(); - - m_h0_RT.resize(m_kk); - m_g0_RT.resize(m_kk); - m_expg0_RT.resize(m_kk); - m_cp0_R.resize(m_kk); - m_s0_R.resize(m_kk); - m_pe.resize(m_kk, 0.0); - m_pp.resize(m_kk); - m_speciesMolarVolume.resize(m_kk); -} - void IdealSolidSolnPhase::setToEquilState(const doublereal* lambda_RT) { const vector_fp& grt = gibbs_RT_ref(); diff --git a/src/thermo/IdealSolnGasVPSS.cpp b/src/thermo/IdealSolnGasVPSS.cpp index ee636f28a..f85245687 100644 --- a/src/thermo/IdealSolnGasVPSS.cpp +++ b/src/thermo/IdealSolnGasVPSS.cpp @@ -233,12 +233,6 @@ void IdealSolnGasVPSS::getPartialMolarVolumes(doublereal* vbar) const getStandardVolumes(vbar); } -void IdealSolnGasVPSS::initThermo() -{ - initLengths(); - VPStandardStateTP::initThermo(); -} - void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT) { updateStandardStateThermo(); @@ -268,15 +262,17 @@ void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT) setState_PX(pres, &m_pp[0]); } -void IdealSolnGasVPSS::initLengths() +bool IdealSolnGasVPSS::addSpecies(shared_ptr spec) { - m_pp.resize(m_kk, 0.0); + bool added = VPStandardStateTP::addSpecies(spec); + if (added) { + m_pp.push_back(0.0); + } + return added; } void IdealSolnGasVPSS::initThermoXML(XML_Node& phaseNode, const std::string& id_) { - IdealSolnGasVPSS::initLengths(); - if (phaseNode.hasChild("thermo")) { XML_Node& thermoNode = phaseNode.child("thermo"); std::string model = thermoNode["model"]; diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index e3e26889d..dd1974a22 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -277,16 +277,20 @@ const vector_fp& LatticePhase::cp_R_ref() const return m_cp0_R; } -void LatticePhase::initThermo() +bool LatticePhase::addSpecies(shared_ptr spec) { - m_Pref = refPressure(); - m_h0_RT.resize(m_kk); - m_g0_RT.resize(m_kk); - m_cp0_R.resize(m_kk); - m_s0_R.resize(m_kk); - m_speciesMolarVolume.resize(m_kk, 0.0); - - ThermoPhase::initThermo(); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + if (m_kk == 1) { + m_Pref = refPressure(); + } + m_h0_RT.push_back(0.0); + m_g0_RT.push_back(0.0); + m_cp0_R.push_back(0.0); + m_s0_R.push_back(0.0); + m_speciesMolarVolume.push_back(0.0); + } + return added; } void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_) diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index 3b8080350..d91e86709 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -360,7 +360,7 @@ void LatticeSolidPhase::installSlavePhases(XML_Node* phaseNode) void LatticeSolidPhase::initThermo() { - initLengths(); + lkstart_.resize(m_lattice.size() + 1); size_t loc = 0; for (size_t n = 0; n < m_lattice.size(); n++) { size_t nsp = m_lattice[n]->nSpecies(); @@ -375,12 +375,14 @@ void LatticeSolidPhase::initThermo() ThermoPhase::initThermo(); } -void LatticeSolidPhase::initLengths() +bool LatticeSolidPhase::addSpecies(shared_ptr spec) { - theta_.resize(m_lattice.size(), 0); - lkstart_.resize(m_lattice.size() + 1); - m_x.resize(m_kk, 0.0); - tmpV_.resize(m_kk, 0.0); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + m_x.push_back(0.0); + tmpV_.push_back(0.0); + } + return added; } void LatticeSolidPhase::_updateThermo() const diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index e1110001d..f1510a57d 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -249,20 +249,21 @@ void MixtureFugacityTP::setStateFromXML(const XML_Node& state) } } -void MixtureFugacityTP::initThermo() +bool MixtureFugacityTP::addSpecies(shared_ptr spec) { - initLengths(); - ThermoPhase::initThermo(); -} - -void MixtureFugacityTP::initLengths() -{ - moleFractions_.resize(m_kk, 0.0); - moleFractions_[0] = 1.0; - m_h0_RT.resize(m_kk, 0.0); - m_cp0_R.resize(m_kk, 0.0); - m_g0_RT.resize(m_kk, 0.0); - m_s0_R.resize(m_kk, 0.0); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + if (m_kk == 1) { + moleFractions_.push_back(1.0); + } else { + moleFractions_.push_back(0.0); + } + m_h0_RT.push_back(0.0); + m_cp0_R.push_back(0.0); + m_g0_RT.push_back(0.0); + m_s0_R.push_back(0.0); + } + return added; } void MixtureFugacityTP::setTemperature(const doublereal temp) @@ -409,12 +410,6 @@ void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublerea setState_TP(t,p); } -void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) -{ - MixtureFugacityTP::initLengths(); - ThermoPhase::initThermoXML(phaseNode, id_); -} - doublereal MixtureFugacityTP::z() const { return pressure() * meanMolecularWeight() / (density() * RT()); diff --git a/src/thermo/MolalityVPSSTP.cpp b/src/thermo/MolalityVPSSTP.cpp index 166fdb0c6..cbadc75d8 100644 --- a/src/thermo/MolalityVPSSTP.cpp +++ b/src/thermo/MolalityVPSSTP.cpp @@ -325,12 +325,8 @@ void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, const std::string& void MolalityVPSSTP::initThermo() { - initLengths(); VPStandardStateTP::initThermo(); - // The solvent defaults to species 0 - setSolvent(0); - // Find the Cl- species m_indexCLM = findCLMIndex(); } @@ -401,20 +397,17 @@ size_t MolalityVPSSTP::findCLMIndex() const return indexCLM; } -// Initialize lengths of local variables after all species have -// been identified. -void MolalityVPSSTP::initLengths() +bool MolalityVPSSTP::addSpecies(shared_ptr spec) { - m_molalities.resize(m_kk); -} - -void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) -{ - initLengths(); - - // The solvent defaults to species 0 - setSolvent(0); - VPStandardStateTP::initThermoXML(phaseNode, id_); + bool added = VPStandardStateTP::addSpecies(spec); + if (added) { + if (m_kk == 1) { + // The solvent defaults to species 0 + setSolvent(0); + } + m_molalities.push_back(0.0); + } + return added; } std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp index ded087320..de071e5d4 100644 --- a/src/thermo/Phase.cpp +++ b/src/thermo/Phase.cpp @@ -834,6 +834,7 @@ bool Phase::addSpecies(shared_ptr spec) { m_y.push_back(0.0); m_ym.push_back(0.0); } + invalidateCache(); return true; } diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 56b03bec3..3c3ee3107 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -34,6 +34,7 @@ RedlichKwongMFTP::RedlichKwongMFTP() : m_b_current(0.0), m_a_current(0.0), NSolns_(0), + Vroot_{0.0, 0.0, 0.0}, dpdV_(0.0), dpdT_(0.0) { @@ -45,6 +46,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string& m_b_current(0.0), m_a_current(0.0), NSolns_(0), + Vroot_{0.0, 0.0, 0.0}, dpdV_(0.0), dpdT_(0.0) { @@ -57,6 +59,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id m_b_current(0.0), m_a_current(0.0), NSolns_(0), + Vroot_{0.0, 0.0, 0.0}, dpdV_(0.0), dpdT_(0.0) { @@ -541,16 +544,6 @@ doublereal RedlichKwongMFTP::critDensity() const return mmw / vc; } -void RedlichKwongMFTP::initThermo() -{ - Vroot_[0] = 0.0; - Vroot_[1] = 0.0; - Vroot_[2] = 0.0; - - initLengths(); - MixtureFugacityTP::initThermo(); -} - void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT) { double tmp, tmp2; @@ -581,27 +574,29 @@ void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT) setState_PX(pres, &m_pp[0]); } -void RedlichKwongMFTP::initLengths() +bool RedlichKwongMFTP::addSpecies(shared_ptr spec) { - a_vec_Curr_.resize(m_kk * m_kk, 0.0); - b_vec_Curr_.resize(m_kk, 0.0); + bool added = MixtureFugacityTP::addSpecies(spec); + if (added) { + a_vec_Curr_.resize(m_kk * m_kk, 0.0); + b_vec_Curr_.push_back(0.0); - a_coeff_vec.resize(2, m_kk * m_kk, 0.0); + a_coeff_vec.resize(2, m_kk * m_kk, 0.0); - m_pc_Species.resize(m_kk, 0.0); - m_tc_Species.resize(m_kk, 0.0); - m_vc_Species.resize(m_kk, 0.0); + m_pc_Species.push_back(0.0); + m_tc_Species.push_back(0.0); + m_vc_Species.push_back(0.0); - m_pp.resize(m_kk, 0.0); - m_tmpV.resize(m_kk, 0.0); - m_partialMolarVolumes.resize(m_kk, 0.0); - dpdni_.resize(m_kk, 0.0); + m_pp.push_back(0.0); + m_tmpV.push_back(0.0); + m_partialMolarVolumes.push_back(0.0); + dpdni_.push_back(0.0); + } + return added; } void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id) { - RedlichKwongMFTP::initLengths(); - // Check the model parameter for the Redlich-Kwong equation of state // two are allowed // RedlichKwong mixture of species, each of which are RK fluids diff --git a/src/thermo/SingleSpeciesTP.cpp b/src/thermo/SingleSpeciesTP.cpp index ccb7ebfe7..a57f476bb 100644 --- a/src/thermo/SingleSpeciesTP.cpp +++ b/src/thermo/SingleSpeciesTP.cpp @@ -20,7 +20,10 @@ namespace Cantera { SingleSpeciesTP::SingleSpeciesTP() : m_press(OneAtm), - m_p0(OneAtm) + m_p0(OneAtm), + m_h0_RT(1), + m_cp0_R(1), + m_s0_R(1) { } @@ -276,25 +279,18 @@ void SingleSpeciesTP::setState_SV(doublereal s, doublereal v, throw CanteraError("setState_SV","no convergence. dt = {}", dt); } -void SingleSpeciesTP::initThermo() +bool SingleSpeciesTP::addSpecies(shared_ptr spec) { - // Make sure there is one and only one species in this phase. - if (nSpecies() != 1) { - throw CanteraError("initThermo", - "stoichiometric substances may only contain one species."); + if (m_kk != 0) { + throw CanteraError("SingleSpeciesTP::addSpecies", + "Stoichiometric substances may only contain one species."); } - - // Resize temporary arrays. - m_h0_RT.resize(1); - m_cp0_R.resize(1); - m_s0_R.resize(1); - - // Make sure the species mole fraction is equal to 1.0; - double x = 1.0; - ThermoPhase::setMoleFractions(&x); - - // Call the base class initThermo object. - ThermoPhase::initThermo(); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + double x = 1.0; + ThermoPhase::setMoleFractions(&x); + } + return added; } void SingleSpeciesTP::_updateThermo() const diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp index c5f9f8794..bcdfab9e2 100644 --- a/src/thermo/SurfPhase.cpp +++ b/src/thermo/SurfPhase.cpp @@ -232,21 +232,22 @@ void SurfPhase::getCp_R_ref(doublereal* cprt) const getCp_R(cprt); } -void SurfPhase::initThermo() +bool SurfPhase::addSpecies(shared_ptr spec) { - ThermoPhase::initThermo(); - m_h0.resize(m_kk); - m_s0.resize(m_kk); - m_cp0.resize(m_kk); - m_mu0.resize(m_kk); - m_work.resize(m_kk); - vector_fp cov(m_kk, 0.0); - cov[0] = 1.0; - setCoverages(cov.data()); - m_logsize.resize(m_kk); - for (size_t k = 0; k < m_kk; k++) { - m_logsize[k] = log(size(k)); + bool added = ThermoPhase::addSpecies(spec); + if (added) { + m_h0.push_back(0.0); + m_s0.push_back(0.0); + m_cp0.push_back(0.0); + m_mu0.push_back(0.0); + m_work.push_back(0.0); + m_logsize.push_back(log(size(m_kk-1))); + if (m_kk == 1) { + vector_fp cov{1.0}; + setCoverages(cov.data()); + } } + return added; } void SurfPhase::setSiteDensity(doublereal n0) diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index 6ce1f801a..a812c1125 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -689,17 +689,11 @@ void ThermoPhase::getReferenceComposition(doublereal* const x) const void ThermoPhase::initThermo() { - // Check to see that there is at least one species defined in the phase - if (m_kk == 0) { - throw CanteraError("ThermoPhase::initThermo()", - "Number of species is equal to zero"); - } // Check to see that all of the species thermo objects have been initialized if (!m_spthermo->ready(m_kk)) { throw CanteraError("ThermoPhase::initThermo()", "Missing species thermo data"); } - xMol_Ref.resize(m_kk, 0.0); } void ThermoPhase::installSlavePhases(XML_Node* phaseNode) { @@ -711,6 +705,7 @@ bool ThermoPhase::addSpecies(shared_ptr spec) if (added) { spec->thermo->validate(spec->name); m_spthermo->install_STIT(m_kk-1, spec->thermo); + xMol_Ref.push_back(0.0); } return added; }