From e7e02ffc5e709a7fab44e2a1be1233007b6070d4 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 17 Dec 2008 17:01:29 +0000 Subject: [PATCH] Changes to exit() and additions of cstdlib --- Cantera/src/thermo/DebyeHuckel.cpp | 25 ++-- Cantera/src/thermo/HMWSoln.cpp | 41 +++---- Cantera/src/thermo/HMWSoln_input.cpp | 6 +- Cantera/src/thermo/MolalityVPSSTP.cpp | 118 +++++++++++++++++- Cantera/src/thermo/MolalityVPSSTP.h | 139 +++++++++++++++++----- Cantera/src/thermo/WaterPropsIAPWS.cpp | 6 +- Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 4 +- Cantera/src/thermo/phasereport.cpp | 2 +- 8 files changed, 268 insertions(+), 73 deletions(-) diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index 20c4292d3..ca4985702 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -26,6 +26,7 @@ #include "WaterProps.h" #include "PDSS_Water.h" #include +#include using namespace std; @@ -997,11 +998,11 @@ namespace Cantera { if (formString == "unity") { m_formGC = 0; printf("exit standardConc = unity not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "molar_volume") { m_formGC = 1; printf("exit standardConc = molar_volume not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "solvent_volume") { m_formGC = 2; } else { @@ -1113,11 +1114,11 @@ namespace Cantera { if (formString == "unity") { m_formGC = 0; printf("exit standardConc = unity not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "molar_volume") { m_formGC = 1; printf("exit standardConc = molar_volume not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "solvent_volume") { m_formGC = 2; } else { @@ -1682,7 +1683,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - exit(-1); + exit(EXIT_FAILURE); } return A; } @@ -1716,7 +1717,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - exit(-1); + exit(EXIT_FAILURE); } return dAdT; } @@ -1750,7 +1751,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - exit(-1); + exit(EXIT_FAILURE); } return d2AdT2; } @@ -1784,7 +1785,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - exit(-1); + exit(EXIT_FAILURE); } return dAdP; } @@ -2176,7 +2177,7 @@ namespace Cantera { default: printf("ERROR\n"); - exit(-1); + exit(EXIT_FAILURE); } /* * Above, we calculated the ln(activitySolvent). Translate that @@ -2321,7 +2322,7 @@ namespace Cantera { default: printf("ERROR\n"); - exit(-1); + exit(EXIT_FAILURE); break; } @@ -2458,7 +2459,7 @@ namespace Cantera { default: printf("ERROR\n"); - exit(-1); + exit(EXIT_FAILURE); break; } } @@ -2596,7 +2597,7 @@ namespace Cantera { default: printf("ERROR\n"); - exit(-1); + exit(EXIT_FAILURE); break; } } diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index 25a6c27a7..75f611cba 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -29,7 +29,8 @@ #include "ThermoFactory.h" #include "WaterProps.h" #include "PDSS_Water.h" -#include +#include +#include namespace Cantera { @@ -389,7 +390,7 @@ namespace Cantera { { if (testProb != 1) { printf("unknown test problem\n"); - std::exit(-1); + exit(EXIT_FAILURE); } constructPhaseFile("HMW_NaCl.xml", ""); @@ -1324,7 +1325,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - std::exit(-1); + exit(EXIT_FAILURE); } return A; } @@ -1359,7 +1360,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - std::exit(-1); + exit(EXIT_FAILURE); } return dAdT; } @@ -1393,7 +1394,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - std::exit(-1); + exit(EXIT_FAILURE); } return dAdP; } @@ -1500,7 +1501,7 @@ namespace Cantera { break; default: printf("shouldn't be here\n"); - std::exit(-1); + exit(EXIT_FAILURE); } return d2AdT2; } @@ -2175,7 +2176,7 @@ namespace Cantera { */ if (m_indexSolvent != 0) { printf("Wrong index solvent value!\n"); - std::exit(-1); + exit(EXIT_FAILURE); } #ifdef DEBUG_MODE @@ -2979,7 +2980,7 @@ namespace Cantera { if (j == (m_kk-1)) { // we should never reach this step printf("logic error 1 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] > 0.0) { /* @@ -3009,7 +3010,7 @@ namespace Cantera { if (j == m_kk-1) { // we should never reach this step printf("logic error 2 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] < 0) { /* @@ -3149,7 +3150,7 @@ namespace Cantera { #endif if (m_indexSolvent != 0) { printf("Wrong index solvent value!\n"); - std::exit(-1); + exit(EXIT_FAILURE); } double d_wateract_dT; @@ -3787,7 +3788,7 @@ namespace Cantera { if (j == (m_kk-1)) { // we should never reach this step printf("logic error 1 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] > 0.0) { /* @@ -3817,7 +3818,7 @@ namespace Cantera { if (j == m_kk-1) { // we should never reach this step printf("logic error 2 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] < 0) { /* @@ -3938,7 +3939,7 @@ namespace Cantera { #endif if (m_indexSolvent != 0) { printf("Wrong index solvent value!\n"); - std::exit(-1); + exit(EXIT_FAILURE); } std::string sni, snj, snk; @@ -4583,7 +4584,7 @@ namespace Cantera { if (j == (m_kk-1)) { // we should never reach this step printf("logic error 1 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] > 0.0) { /* @@ -4613,7 +4614,7 @@ namespace Cantera { if (j == m_kk-1) { // we should never reach this step printf("logic error 2 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] < 0) { /* @@ -4754,7 +4755,7 @@ namespace Cantera { #endif if (m_indexSolvent != 0) { printf("Wrong index solvent value!\n"); - std::exit(-1); + exit(EXIT_FAILURE); } double d_wateract_dP; @@ -5397,7 +5398,7 @@ namespace Cantera { if (j == (m_kk-1)) { // we should never reach this step printf("logic error 1 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] > 0.0) { /* @@ -5428,7 +5429,7 @@ namespace Cantera { if (j == m_kk-1) { // we should never reach this step printf("logic error 2 in Step 9 of hmw_act"); - std::exit(1); + exit(EXIT_FAILURE); } if (charge[k] < 0) { /* @@ -5614,13 +5615,13 @@ namespace Cantera { #ifdef DEBUG_MODE if (i > 4 || j > 4) { printf("we shouldn't be here\n"); - std::exit(-1); + exit(EXIT_FAILURE); } #endif if ((i == 0) || (j == 0)) { printf("ERROR calc_thetas called with one species being neutral\n"); - std::exit(-1); + exit(EXIT_FAILURE); } /* diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index 691a836dc..068fb9968 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -21,10 +21,10 @@ #include "WaterProps.h" #include "PDSS_Water.h" #include +#include using namespace std; - namespace Cantera { @@ -1011,11 +1011,11 @@ namespace Cantera { if (formString == "unity") { m_formGC = 0; printf("exit standardConc = unity not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "molar_volume") { m_formGC = 1; printf("exit standardConc = molar_volume not done\n"); - exit(-1); + exit(EXIT_FAILURE); } else if (formString == "solvent_volume") { m_formGC = 2; } else { diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index ab6b06274..4926a5bf6 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -41,6 +41,8 @@ namespace Cantera { MolalityVPSSTP::MolalityVPSSTP() : VPStandardStateTP(), m_indexSolvent(0), + m_pHScalingType(PHSCALE_PITZER), + m_indexCLM(-1), m_weightSolvent(18.01528), m_xmolSolventMIN(0.01), m_Mnaught(18.01528E-3) @@ -62,6 +64,8 @@ namespace Cantera { MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) : VPStandardStateTP(), m_indexSolvent(b.m_indexSolvent), + m_pHScalingType(b.m_pHScalingType), + m_indexCLM(b.m_indexCLM), m_xmolSolventMIN(b.m_xmolSolventMIN), m_Mnaught(b.m_Mnaught), m_molalities(b.m_molalities) @@ -80,6 +84,8 @@ namespace Cantera { if (&b != this) { VPStandardStateTP::operator=(b); m_indexSolvent = b.m_indexSolvent; + m_pHScalingType = b.m_pHScalingType; + m_indexCLM = b.m_indexCLM; m_weightSolvent = b.m_weightSolvent; m_xmolSolventMIN = b.m_xmolSolventMIN; m_Mnaught = b.m_Mnaught; @@ -391,6 +397,24 @@ namespace Cantera { return cAC_CONVENTION_MOLALITY; } + void MolalityVPSSTP::getActivityConcentrations(doublereal* c) const { + err("getActivityConcentrations"); + } + + doublereal MolalityVPSSTP::standardConcentration(int k) const { + err("standardConcentration"); + return -1.0; + } + + doublereal MolalityVPSSTP::logStandardConc(int k) const { + err("logStandardConc"); + return -1.0; + } + + void MolalityVPSSTP::getActivities(doublereal* ac) const { + err("getActivities"); + } + /* * Get the array of non-dimensional activity coefficients at * the current solution temperature, pressure, and @@ -416,7 +440,12 @@ namespace Cantera { } } } - + + + void MolalityVPSSTP::getMolalityActivityCoefficients(doublereal *acMolality) const { + err("getMolalityActivityCoefficients"); + } + /* * osmotic coefficient: * @@ -452,6 +481,15 @@ namespace Cantera { return oc; } + + void MolalityVPSSTP::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + /* * ------------ Partial Molar Properties of the Solution ------------ */ @@ -495,7 +533,12 @@ namespace Cantera { if (i == 5) uA[5] = 0.0; } } - + + void MolalityVPSSTP::setToEquilState(const doublereal* lambda_RT) { + updateStandardStateThermo(); + err("setToEquilState"); + } + /* * Set the thermodynamic state. */ @@ -561,6 +604,77 @@ namespace Cantera { setSolvent(0); } + // Returns the index of the Cl- species. + /* + * The Cl- species is special in the sense that it's single ion + * molalality-based activity coefficient is used in the specification + * of the pH scale for single ions. Therefore, we need to know + * what species index is Cl-. If the species isn't in the species + * list then this routine returns -1, and we can't use the NBS + * pH scale. + * + * Right now we use a restrictive interpretation. The species + * must be named "Cl-". It must consist of exactly one Cl and one E + * atom. + */ + int MolalityVPSSTP::findCLMIndex() const { + int indexCLM = -1; + int eCl = -1; + int eE = -1; + int ne= nElements(); + string sn; + for (int e = 0; e < ne; e++) { + sn = elementName(e); + if (sn == "Cl" || sn == "CL") { + eCl = e; + break; + } + } + // We have failed if we can't find the Cl element index + if (eCl == -1) { + return -1; + } + for (int e = 0; e < ne; e++) { + sn = elementName(e); + if (sn == "E" || sn == "e") { + eE = e; + break; + } + } + // We have failed if we can't find the E element index + if (eE == -1) { + return -1; + } + for (int k = 1; k < m_kk; k++) { + doublereal nCl = nAtoms(k, eCl); + if (nCl != 1.0) { + continue; + } + doublereal nE = nAtoms(k, eE); + if (nE != -1.0) { + continue; + } + for (int e = 0; e < ne; e++) { + if (e != eE && e != eCl) { + doublereal nA = nAtoms(k, e); + if (nA != 0.0) { + continue; + } + } + } + sn = speciesName(k); + if (sn != "Cl-" && sn != "CL-") { + continue; + } + + indexCLM = k; + break; + } + return indexCLM; + } + + // Initialize lengths of local variables after all species have + // been identified. void MolalityVPSSTP::initLengths() { m_kk = nSpecies(); m_molalities.resize(m_kk); diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index a6f1d1b3b..d36e74bb9 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -138,8 +138,8 @@ namespace Cantera { * return the molality-based quantities. Also all functions which return * activities return the molality-based activities. The reason for this convention * has been discussed in supporting memos. However, it's important because the - * term in the equation above is non-trivial. For example it's equal to 2.38 kcal gmol-1 - * for water at 298 K. + * term in the equation above is non-trivial. For example it's equal + * to 2.38 kcal gmol-1 for water at 298 K. * * * In order to prevent a singularity, this class includes the concept of a minimum @@ -222,6 +222,14 @@ namespace Cantera { virtual int eosType() const { return 0; } + //! Set the pH scale, which determines the scale for single-ion activity + //! coefficients. + /*! + * Single ion activity coefficients are not unique in terms of the + * representing actual measureable quantities. + */ + void setpHScale(const int pHscaleType); + /** * @} * @name Molar Thermodynamic Properties @@ -417,9 +425,7 @@ namespace Cantera { * units depend upon the implementation of the * reaction rate expressions within the phase. */ - virtual void getActivityConcentrations(doublereal* c) const { - err("getActivityConcentrations"); - } + virtual void getActivityConcentrations(doublereal* c) const; /** * The standard concentration \f$ C^0_k \f$ used to normalize @@ -434,10 +440,7 @@ namespace Cantera { * * @param k species index. Defaults to zero. */ - virtual doublereal standardConcentration(int k=0) const { - err("standardConcentration"); - return -1.0; - } + virtual doublereal standardConcentration(int k=0) const; /** * Returns the natural logarithm of the standard @@ -445,10 +448,7 @@ namespace Cantera { * * @param k species index */ - virtual doublereal logStandardConc(int k=0) const { - err("logStandardConc"); - return -1.0; - } + virtual doublereal logStandardConc(int k=0) const; /** * Returns the units of the standard and generalized @@ -492,9 +492,7 @@ namespace Cantera { * * @param ac Output vector of molality-based activities. Length: m_kk. */ - virtual void getActivities(doublereal* ac) const { - err("getActivities"); - } + virtual void getActivities(doublereal* ac) const; //! Get the array of non-dimensional activity coefficients at //! the current solution temperature, pressure, and solution concentration. @@ -540,9 +538,7 @@ namespace Cantera { * @param acMolality Output vector containing the molality based activity coefficients. * length: m_kk. */ - virtual void getMolalityActivityCoefficients(doublereal *acMolality) const { - err("getMolalityActivityCoefficients"); - } + virtual void getMolalityActivityCoefficients(doublereal *acMolality) const; //! Calculate the osmotic coefficient /*! @@ -576,13 +572,7 @@ namespace Cantera { * @param mu output vector containing the species electrochemical potentials. * Length: m_kk. */ - void getElectrochemPotentials(doublereal* mu) const { - getChemPotentials(mu); - double ve = Faraday * electricPotential(); - for (int k = 0; k < m_kk; k++) { - mu[k] += ve*charge(k); - } - } + void getElectrochemPotentials(doublereal* mu) const; //@} @@ -640,10 +630,7 @@ namespace Cantera { * @param lambda_RT Input vector containing the dimensionless * element potentials. */ - virtual void setToEquilState(const doublereal* lambda_RT) { - updateStandardStateThermo(); - err("setToEquilState"); - } + virtual void setToEquilState(const doublereal* lambda_RT); //@} @@ -753,6 +740,23 @@ namespace Cantera { virtual std::string report(bool show_thermo = true) const; private: + //! Returns the index of the Cl- species. + /*! + * The Cl- species is special in the sense that it's single ion + * molalality-based activity coefficient is used in the specification + * of the pH scale for single ions. Therefore, we need to know + * what species index is Cl-. If the species isn't in the species + * list then this routine returns -1, and we can't use the NBS + * pH scale. + * + * Right now we use a restrictive interpretation. The species + * must be named "Cl-". It must consist of exactly one Cl and one E + * atom. + */ + virtual int findCLMIndex() const; + + //! Initialize lengths of local variables after all species have + //! been identified. void initLengths(); protected: @@ -763,6 +767,24 @@ namespace Cantera { */ int m_indexSolvent; + //! Scaling to be used for output of single-ion species activity + //! coefficients. + /*! + * Index of the species to be used in the single-ion scaling + * law. This is the indentity of the Cl- species for the PHSCALE_NBS + * scaling. + * Either PHSCALE_PITZER or PHSCALE_NBS + */ + int m_pHScalingType; + + //! Index of the phScale species + /*! + * Index of the species to be used in the single-ion scaling + * law. This is the indentity of the Cl- species for the PHSCALE_NBS + * scaling + */ + int m_indexCLM; + //! Molecular weight of the Solvent doublereal m_weightSolvent; @@ -790,11 +812,68 @@ namespace Cantera { * units are (kg/kmol) */ mutable vector_fp m_molalities; + private: + //! Error function + /*! + * Print an error string and exit + * + * @param msg Message to be printed + */ doublereal err(std::string msg) const; }; + + + //! Scale to be used for the output of single-ion activity coefficients + //! is that used by Pitzer. + /*! + * This is the internal scale used within the code. One property is that + * the activity coefficients for the cation and anion of a single salt + * will be equal. This scale is the one presumed by the formulation of the + * single-ion activity coefficients described in this report. + * + * Activity coefficients for species k may be altered between scales s1 to s2 + * using the following formula + * + * \f[ + * ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) + * + \frac{z_k}{z_j} \left( ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right) + * \f] + * + * where j is any one species. + * + * + */ + const int PHSCALE_PITZER = 0; + + //! Scale to be used for evaluation of single-ion activity coefficients + //! is that used by the NBS standard for evaluation of the pH variable. + /*! + * This is not the internal scale used within the code. + * + * Activity coefficients for species k may be altered between scales s1 to s2 + * using the following formula + * + * \f[ + * ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) + * + \frac{z_k}{z_j} \left( ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right) + * \f] + * + * where j is any one species. For the NBS scale, j is equal to the Cl- species + * and + * + * \f[ + * ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}} + * \f] + * + * This is the NBS pH scale, which is used in all conventional pH + * measurements. and is based on the Bates-Guggenheim quations. + * + */ + const int PHSCALE_NBS = 1; + } #endif diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index a39c4d246..fcb07119a 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -16,9 +16,9 @@ #include "WaterPropsIAPWS.h" #include "ctexceptions.h" #include "stringUtils.h" -#include -#include -#include +#include +#include +#include /* * Critical Point values of water in mks units */ diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index 548bbd139..5dde4a863 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -14,8 +14,8 @@ #include "WaterPropsIAPWSphi.h" -#include -#include +#include +#include /* * Critical Point values in mks units: Note, these aren't used in this diff --git a/Cantera/src/thermo/phasereport.cpp b/Cantera/src/thermo/phasereport.cpp index 9b9986245..a5be76217 100644 --- a/Cantera/src/thermo/phasereport.cpp +++ b/Cantera/src/thermo/phasereport.cpp @@ -19,7 +19,7 @@ #include "ThermoPhase.h" #include "PureFluidPhase.h" -#include +#include #include "mix_defs.h" using namespace std;