From e41f43300f7703b26938724e5f1cccf21a8aa2b9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 9 Mar 2012 22:55:34 +0000 Subject: [PATCH] Folded class State into class Phase --- include/cantera/thermo/DebyeHuckel.h | 2 +- include/cantera/thermo/HMWSoln.h | 2 +- include/cantera/thermo/Phase.h | 386 +++++++++++++++++-- include/cantera/thermo/RedlichKwongMFTP.h | 2 +- include/cantera/thermo/State.h | 449 ---------------------- include/cantera/thermo/SurfPhase.h | 4 +- src/thermo/DebyeHuckel.cpp | 4 +- src/thermo/GibbsExcessVPSSTP.cpp | 14 +- src/thermo/HMWSoln.cpp | 6 +- src/thermo/IdealMolalSoln.cpp | 6 +- src/thermo/IdealSolidSolnPhase.cpp | 14 +- src/thermo/IdealSolnGasVPSS.cpp | 6 +- src/thermo/IonsFromNeutralVPSSTP.cpp | 2 +- src/thermo/LatticePhase.cpp | 14 +- src/thermo/LatticeSolidPhase.cpp | 6 +- src/thermo/MixtureFugacityTP.cpp | 34 +- src/thermo/Phase.cpp | 247 ++++++++++-- src/thermo/RedlichKwongMFTP.cpp | 6 +- src/thermo/State.cpp | 311 --------------- src/thermo/SurfPhase.cpp | 6 +- src/thermo/VPStandardStateTP.cpp | 2 +- src/thermo/WaterSSTP.cpp | 8 +- 22 files changed, 646 insertions(+), 885 deletions(-) delete mode 100644 include/cantera/thermo/State.h delete mode 100644 src/thermo/State.cpp diff --git a/include/cantera/thermo/DebyeHuckel.h b/include/cantera/thermo/DebyeHuckel.h index e2fc09ed6..c852d159e 100644 --- a/include/cantera/thermo/DebyeHuckel.h +++ b/include/cantera/thermo/DebyeHuckel.h @@ -832,7 +832,7 @@ public: * the value propagates to underlying objects, such as * the water standard state model. * - * @todo Make State::setTemperature a virtual function + * @todo Make Phase::setTemperature a virtual function * * @param temp Temperature in kelvin */ diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h index f2af6f2c9..e84a1b32d 100644 --- a/include/cantera/thermo/HMWSoln.h +++ b/include/cantera/thermo/HMWSoln.h @@ -1533,7 +1533,7 @@ public: * the value propagates to underlying objects, such as * the water standard state model. * - * @todo Make State::setTemperature a virtual function + * @todo Make Phase::setTemperature a virtual function * * @param temp Temperature in kelvin */ diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 10c74cca3..a7b1a4866 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -1,17 +1,16 @@ /** * @file Phase.h * - * Header file for class, Phase, which contains functions for - * setting the state of a phase, and for referencing species by - * name, and also contains text for the module phases (see \ref - * phases and class \link Cantera::Phase Phase\endlink). + * Header file for class, Phase, which manages the independent variables + * of temperature, mass density, and species mass/mole fraction that define + * the thermodynamic state. Also contains functions for managing the species + * and elements in the phase. (see \ref phases) */ // Copyright 2001 California Institute of Technology #ifndef CT_PHASE_H #define CT_PHASE_H -#include "State.h" #include "Constituents.h" #include "cantera/base/vec_functions.h" @@ -27,6 +26,37 @@ namespace Cantera * classes, it implements methods that allow referencing a species * by name. * + * Manages the independent variables of temperature, mass density, + * and species mass/mole fraction that define the thermodynamic + * state. + * Class State stores just enough information about a + * multicomponent solution to specify its intensive thermodynamic + * state. It stores values for the temperature, mass density, and + * an array of species mass fractions. It also stores an array of + * species molecular weights, which are used to convert between + * mole and mass representations of the composition. These are the + * \e only properties of the species that class State knows about. + * For efficiency in mass/mole conversion, the vector of mass + * fractions divided by molecular weight \f$ Y_k/M_k \f$ is also + * stored. + * + * Class State is not usually used directly in application + * programs. Its primary use is as a base class for class + * Phase. Class State has no virtual methods, and none of its + * methods are meant to be overloaded. However, this is one exception. + * If the phase is incompressible, then the density must be replaced + * by the pressure as the independent variable. In this case, functions + * such as setMassFraction within the class %State must actually now + * calculate the density (at constant T and P) instead of leaving + * it alone as befits an independent variable. Threfore, these type + * of functions are virtual functions and need to be overloaded + * for incompressible phases. Note, for almost incompressible phases + * (or phases which utilize standard states based on a T and P) this + * may be advantageous as well, and they need to overload these functions + * too. + * + * @ingroup phases + * * Class Phase derives from both classes * Constituents and State. In addition to the methods of those two * classes, it implements methods that allow referencing a species @@ -81,11 +111,9 @@ namespace Cantera * * @ingroup phases */ -class Phase : public Constituents, public State +class Phase : public Constituents { - public: - /// Default constructor. Phase(); @@ -258,8 +286,7 @@ public: * * @param t Temperature in kelvin * @param dens Density (kg/m^3) - * @param x vector of species mole fractions. - * Length is equal to m_kk + * @param x vector of species mole fractions, length m_kk */ void setState_TRX(doublereal t, doublereal dens, const doublereal* x); @@ -282,8 +309,7 @@ public: * * @param t Temperature in kelvin * @param dens Density (kg/m^3) - * @param y vector of species mass fractions. - * Length is equal to m_kk + * @param y vector of species mass fractions, length m_kk */ void setState_TRY(doublereal t, doublereal dens, const doublereal* y); @@ -305,8 +331,7 @@ public: * * @param t Temperature in kelvin * @param n molar density (kmol/m^3) - * @param x vector of species mole fractions. - * Length is equal to m_kk + * @param x vector of species mole fractions, length m_kk */ void setState_TNX(doublereal t, doublereal n, const doublereal* x); @@ -320,32 +345,28 @@ public: //! Set the internally stored temperature (K) and mole fractions. /*! * @param t Temperature in kelvin - * @param x vector of species mole fractions. - * Length is equal to m_kk + * @param x vector of species mole fractions, length m_kk */ void setState_TX(doublereal t, doublereal* x); //! Set the internally stored temperature (K) and mass fractions. /*! * @param t Temperature in kelvin - * @param y vector of species mass fractions. - * Length is equal to m_kk + * @param y vector of species mass fractions, length m_kk */ void setState_TY(doublereal t, doublereal* y); //! Set the density (kg/m^3) and mole fractions. /*! * @param rho Density (kg/m^3) - * @param x vector of species mole fractions. - * Length is equal to m_kk + * @param x vector of species mole fractions, length m_kk */ void setState_RX(doublereal rho, doublereal* x); //! Set the density (kg/m^3) and mass fractions. /*! * @param rho Density (kg/m^3) - * @param y vector of species mass fractions. - * Length is equal to m_kk + * @param y vector of species mass fractions, length m_kk */ void setState_RY(doublereal rho, doublereal* y); @@ -374,8 +395,7 @@ public: void getMolecularWeights(doublereal* weights) const; /** - * Return a const reference to the internal vector of - * molecular weights. + * Return a const reference to the internal vector of molecular weights. */ const vector_fp& molecularWeights() const; @@ -389,8 +409,7 @@ public: //! Return the mole fraction of a single species /*! - * @param k String name of the species - * + * @param k species index * @return Mole fraction of the species */ doublereal moleFraction(size_t k) const; @@ -398,27 +417,141 @@ public: //! Return the mole fraction of a single species /*! * @param name String name of the species - * * @return Mole fraction of the species */ doublereal moleFraction(std::string name) const; //! Return the mass fraction of a single species /*! - * @param k String name of the species - * - * @return Mass Fraction of the species + * @param k species index + * @return Mass fraction of the species */ doublereal massFraction(size_t k) const; //! Return the mass fraction of a single species /*! * @param name String name of the species - * * @return Mass Fraction of the species */ doublereal massFraction(std::string name) const; + //@} + /// @name Composition + //@{ + + //! Get the species mole fraction vector. + /*! + * @param x On return, x contains the mole fractions. Must have a + * length greater than or equal to the number of species. + */ + void getMoleFractions(doublereal* const x) const; + + //! Set the mole fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param x Array of unnormalized mole fraction values (input). + * Must have a length greater than or equal to the number of + * species, m_kk. There is no restriction + * on the sum of the mole fraction vector. Internally, + * the State object will normalize this vector before + * storing its contents. + */ + virtual void setMoleFractions(const doublereal* const x); + + /** + * Set the mole fractions to the specified values without + * normalizing. This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of + * equations. + * + * @param x Input vector of mole fractions. + * Length is m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal* const x); + + //! Get the species mass fractions. + /*! + * @param y On return, y contains the mass fractions. Array \a y must have a length + * greater than or equal to the number of species. + */ + void getMassFractions(doublereal* const y) const; + + //! Returns a read-only pointer to the start of the massFraction array + /*! + * @return returns a pointer to a vector of doubles of length m_kk. + */ + const doublereal* massFractions() const { + return &m_y[0]; + } + + //! Set the mass fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param y Array of unnormalized mass fraction values (input). + * Must have a length greater than or equal to the number of species. + * Input vector of mass fractions. There is no restriction + * on the sum of the mass fraction vector. Internally, + * the State object will normalize this vector before + * storing its contents. + * Length is m_kk. + */ + virtual void setMassFractions(const doublereal* const y); + + //! Set the mass fractions to the specified values without normalizing. + /*! + * This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of equations. + * + * @param y Input vector of mass fractions. + * Length is m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal* const y); + + /** + * Get the species concentrations (kmol/m^3). @param c On + * return, \a c contains the concentrations for all species. + * Array \a c must have a length greater than or equal to the + * number of species. + */ + void getConcentrations(doublereal* const c) const; + + /** + * Concentration of species k. If k is outside the valid + * range, an exception will be thrown. + * + * @param k Index of species + */ + doublereal concentration(const size_t k) const; + + //! Set the concentrations to the specified values within the + //! phase. + /*! + * We set the concentrations here and therefore we set the + * overall density of the phase. We hold the temperature constant + * during this operation. Therefore, we have possibly changed + * the pressure of the phase by calling this routine. + * + * @param conc The input vector to this routine is in dimensional + * units. For volumetric phases c[k] is the + * concentration of the kth species in kmol/m3. + * For surface phases, c[k] is the concentration + * in kmol/m2. The length of the vector is the number + * of species in the phase. + */ + virtual void setConcentrations(const doublereal* const conc); + + /** + * Returns a read-only pointer to the start of the + * moleFraction/MW array. This array is the array of mole + * fractions, each divided by the mean molecular weight. + */ + const doublereal* moleFractdivMMW() const; + + //@} + + /** * Charge density [C/m^3]. */ @@ -440,6 +573,113 @@ public: m_ndim = ndim; } + /// @name Thermodynamic Properties + /// Class Phase only stores enough thermodynamic data to + /// specify the state. In addition to composition information, + /// it stores the temperature and mass density. + //@{ + + //! Temperature (K). + /*! + * @return Returns the temperature of the phase + */ + doublereal temperature() const { + return m_temp; + } + + //! Density (kg/m^3). + /*! + * @return Returns the density of the phase + */ + virtual doublereal density() const { + return m_dens; + } + + //! Molar density (kmol/m^3). + /*! + * @return Returns the molar density of the phase + */ + doublereal molarDensity() const; + + //! Molar volume (m^3/kmol). + /*! + * @return Returns the molar volume of the phase + */ + doublereal molarVolume() const; + + //! Set the internally stored density (kg/m^3) of the phase + /*! + * Note the density of a phase is an indepedent variable. + * + * @param density Input density (kg/m^3). + */ + virtual void setDensity(const doublereal density) { + m_dens = density; + } + + //! Set the internally stored molar density (kmol/m^3) of the phase. + /*! + * @param molarDensity Input molar density (kmol/m^3). + */ + virtual void setMolarDensity(const doublereal molarDensity); + + //! Set the temperature (K). + /*! + * This function sets the internally stored temperature of the phase. + * + * @param temp Temperature in kelvin + */ + virtual void setTemperature(const doublereal temp) { + m_temp = temp; + } + //@} + + /// @name Mean Properties + //@{ + /** + * Evaluate the mole-fraction-weighted mean of Q: + * \f[ \sum_k X_k Q_k. \f] + * Array Q should contain pure-species molar property + * values. + * + * @param Q input vector of length m_kk that is to be averaged. + * @return + * mole-freaction-weighted mean of Q + */ + doublereal mean_X(const doublereal* const Q) const; + + /** + * Evaluate the mass-fraction-weighted mean of Q: + * \f[ \sum_k Y_k Q_k \f] + * + * @param Q Array Q contains a vector of species property values in mass units. + * @return + * Return value containing the mass-fraction-weighted mean of Q. + */ + doublereal mean_Y(const doublereal* const Q) const; + + /** + * The mean molecular weight. Units: (kg/kmol) + */ + doublereal meanMolecularWeight() const { + return m_mmw; + } + + //! Evaluate \f$ \sum_k X_k \log X_k \f$. + /*! + * @return + * returns the indicated sum. units are dimensionless. + */ + doublereal sum_xlogx() const; + + //! Evaluate \f$ \sum_k X_k \log Q_k \f$. + /*! + * @param Q Vector of length m_kk to take the log average of + * @return Returns the indicated sum. + */ + doublereal sum_xlogQ(doublereal* const Q) const; + //@} + /** * Finished adding species, prepare to use them for calculation * of mixture properties. @@ -448,8 +688,41 @@ public: virtual bool ready() const; + //! Return the State Mole Fraction Number + DEPRECATED(int stateMFNumber() const) { + return m_stateNum; + } + + //! Every time the mole fractions have changed, this routine + //! will increment the stateMFNumber + /*! + * @param forceChange If this is true then the stateMFNumber always + * changes. This defaults to false. + * @deprecated + */ + DEPRECATED(void stateMFChangeCalc(bool forceChange = false)); protected: + /** + * @internal + * Initialize. Make a local copy of the vector of + * molecular weights, and resize the composition arrays to + * the appropriate size. The only information an instance of + * State has about the species is their molecular weights. + * + * @param mw Vector of molecular weights of the species. + */ + void init(const vector_fp& mw); + + //! Set the molecular weight of a single species to a given value + /*! + * @param k id of the species + * @param mw Molecular Weight (kg kmol-1) + */ + void setMolecularWeight(const int k, const double mw) { + m_molwts[k] = mw; + m_rmolwts[k] = 1.0/mw; + } /** * m_kk = Number of species in the phase. @internal m_kk is a @@ -500,6 +773,57 @@ private: * names and unique phases within a Cantera problem. */ std::string m_name; + + /** + * Temperature. This is an independent variable + * units = Kelvin + */ + doublereal m_temp; + + /** + * Density. This is an independent variable except in + * the incompressible degenerate case. Thus, + * the pressure is determined from this variable + * not the other way round. + * units = kg m-3 + */ + doublereal m_dens; + + /** + * m_mmw is the mean molecular weight of the mixture + * (kg kmol-1) + */ + doublereal m_mmw; + + /** + * m_ym[k] = mole fraction of species k divided by the + * mean molecular weight of mixture. + */ + mutable vector_fp m_ym; + + /** + * m_y[k] = mass fraction of species k + */ + mutable vector_fp m_y; + + /** + * m_molwts[k] = molecular weight of species k (kg kmol-1) + */ + vector_fp m_molwts; + + /** + * m_rmolwts[k] = inverse of the molecular weight of species k + * units = kmol kg-1. + */ + vector_fp m_rmolwts; + + //! State Change variable + /*! + * Whenever the mole fraction vector changes, this int is + * incremented. + * @deprecated + */ + int m_stateNum; }; //! typedef for the base Phase class diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h index e16504c22..074869715 100644 --- a/include/cantera/thermo/RedlichKwongMFTP.h +++ b/include/cantera/thermo/RedlichKwongMFTP.h @@ -220,7 +220,7 @@ protected: * function sets the temperature, and makes sure that * the value propagates to underlying objects * - * @todo Make State::setTemperature a virtual function + * @todo Make Phase::setTemperature a virtual function * * @param temp Temperature in kelvin */ diff --git a/include/cantera/thermo/State.h b/include/cantera/thermo/State.h deleted file mode 100644 index fa5f38c3f..000000000 --- a/include/cantera/thermo/State.h +++ /dev/null @@ -1,449 +0,0 @@ -/** - * @file State.h Header for the class State, that manages the - * independent variables of temperature, mass density, and species - * mass/mole fraction that define the thermodynamic state (see \ref - * phases and class \link Cantera::State State\endlink). - */ - -/* - * Copyright 2001-2003 California Institute of Technology - * See file License.txt for licensing information - */ - -#ifndef CT_STATE2_H -#define CT_STATE2_H - -#include "cantera/base/ct_defs.h" -#include "cantera/base/utilities.h" - -namespace Cantera -{ - -//! Manages the independent variables of temperature, mass density, -//! and species mass/mole fraction that define the thermodynamic -//! state. -/*! - * Class State stores just enough information about a - * multicomponent solution to specify its intensive thermodynamic - * state. It stores values for the temperature, mass density, and - * an array of species mass fractions. It also stores an array of - * species molecular weights, which are used to convert between - * mole and mass representations of the composition. These are the - * \e only properties of the species that class State knows about. - * For efficiency in mass/mole conversion, the vector of mass - * fractions divided by molecular weight \f$ Y_k/M_k \f$ is also - * stored. - * - * Class State is not usually used directly in application - * programs. Its primary use is as a base class for class - * Phase. Class State has no virtual methods, and none of its - * methods are meant to be overloaded. However, this is one exception. - * If the phase is incompressible, then the density must be replaced - * by the pressure as the independent variable. In this case, functions - * such as setMassFraction within the class %State must actually now - * calculate the density (at constant T and P) instead of leaving - * it alone as befits an independent variable. Threfore, these type - * of functions are virtual functions and need to be overloaded - * for incompressible phases. Note, for almost incompressible phases - * (or phases which utilize standard states based on a T and P) this - * may be advantageous as well, and they need to overload these functions - * too. - * - * @ingroup phases - */ -class State -{ - -public: - - /** - * Constructor. - */ - State(); - - /** - * Destructor. Since no memory is allocated by methods of this - * class, the destructor does nothing. - */ - virtual ~State(); - - /** - * Copy Constructor for the State Class - * - * @param right Reference to the class to be copied. - */ - State(const State& right); - - /** - * Assignment operator for the state class. - * - * @param right Reference to the class to be copied. - */ - State& operator=(const State& right); - - - /// @name Species Information - /// - /// The only thing class State knows about the species is their - /// molecular weights. - //@{ - - /// Return a read-only reference to the array of molecular - /// weights. - const vector_fp& molecularWeights() const { - return m_molwts; - } - - - //@} - /// @name Composition - //@{ - - - //! Get the species mole fraction vector. - /*! - * @param x On return, x contains the mole fractions. Must have a - * length greater than or equal to the number of species. - */ - void getMoleFractions(doublereal* const x) const; - - - //! The mole fraction of species k. - /*! - * If k is ouside the valid - * range, an exception will be thrown. Note that it is - * somewhat more efficent to call getMoleFractions if the - * mole fractions of all species are desired. - * @param k species index - */ - doublereal moleFraction(const size_t k) const; - - //! Set the mole fractions to the specified values, and then - //! normalize them so that they sum to 1.0. - /*! - * @param x Array of unnormalized mole fraction values (input). - * Must have a length greater than or equal to the number of - * species, m_kk. There is no restriction - * on the sum of the mole fraction vector. Internally, - * the State object will normalize this vector before - * storing its contents. - */ - virtual void setMoleFractions(const doublereal* const x); - - /** - * Set the mole fractions to the specified values without - * normalizing. This is useful when the normalization - * condition is being handled by some other means, for example - * by a constraint equation as part of a larger set of - * equations. - * - * @param x Input vector of mole fractions. - * Length is m_kk. - */ - virtual void setMoleFractions_NoNorm(const doublereal* const x); - - //! Get the species mass fractions. - /*! - * @param y On return, y contains the mass fractions. Array \a y must have a length - * greater than or equal to the number of species. - */ - void getMassFractions(doublereal* const y) const; - - //! Mass fraction of species k. - /*! - * If k is outside the valid range, an exception will be thrown. Note that it is - * somewhat more efficent to call getMassFractions if the mass fractions of all species are desired. - * - * @param k species index - */ - doublereal massFraction(const size_t k) const; - - //! Set the mass fractions to the specified values, and then - //! normalize them so that they sum to 1.0. - /*! - * @param y Array of unnormalized mass fraction values (input). - * Must have a length greater than or equal to the number of species. - * Input vector of mass fractions. There is no restriction - * on the sum of the mass fraction vector. Internally, - * the State object will normalize this vector before - * storing its contents. - * Length is m_kk. - */ - virtual void setMassFractions(const doublereal* const y); - - //! Set the mass fractions to the specified values without normalizing. - /*! - * This is useful when the normalization - * condition is being handled by some other means, for example - * by a constraint equation as part of a larger set of equations. - * - * @param y Input vector of mass fractions. - * Length is m_kk. - */ - virtual void setMassFractions_NoNorm(const doublereal* const y); - - /** - * Get the species concentrations (kmol/m^3). @param c On - * return, \a c contains the concentrations for all species. - * Array \a c must have a length greater than or equal to the - * number of species. - */ - void getConcentrations(doublereal* const c) const; - - /** - * Concentration of species k. If k is outside the valid - * range, an exception will be thrown. - * - * @param k Index of species - */ - doublereal concentration(const size_t k) const; - - //! Set the concentrations to the specified values within the - //! phase. - /*! - * We set the concentrations here and therefore we set the - * overall density of the phase. We hold the temperature constant - * during this operation. Therefore, we have possibly changed - * the pressure of the phase by calling this routine. - * - * @param conc The input vector to this routine is in dimensional - * units. For volumetric phases c[k] is the - * concentration of the kth species in kmol/m3. - * For surface phases, c[k] is the concentration - * in kmol/m2. The length of the vector is the number - * of species in the phase. - */ - virtual void setConcentrations(const doublereal* const conc); - - //! Returns a read-only pointer to the start of the - //! massFraction array - /*! - * The pointer returned is readonly - * @return returns a pointer to a vector of doubles of length m_kk. - */ - const doublereal* massFractions() const { - return &m_y[0]; - } - - /** - * Returns a read-only pointer to the start of the - * moleFraction/MW array. This array is the array of mole - * fractions, each divided by the mean molecular weight. - */ - const doublereal* moleFractdivMMW() const; - - //@} - - /// @name Mean Properties - //@{ - /** - * Evaluate the mole-fraction-weighted mean of Q: - * \f[ \sum_k X_k Q_k. \f] - * Array Q should contain pure-species molar property - * values. - * - * @param Q input vector of length m_kk that is to be averaged. - * @return - * mole-freaction-weighted mean of Q - */ - doublereal mean_X(const doublereal* const Q) const; - - /** - * Evaluate the mass-fraction-weighted mean of Q: - * \f[ \sum_k Y_k Q_k \f] - * - * @param Q Array Q contains a vector of species property values in mass units. - * @return - * Return value containing the mass-fraction-weighted mean of Q. - */ - doublereal mean_Y(const doublereal* const Q) const; - - /** - * The mean molecular weight. Units: (kg/kmol) - */ - doublereal meanMolecularWeight() const { - return m_mmw; - } - - //! Evaluate \f$ \sum_k X_k \log X_k \f$. - /*! - * @return - * returns the indicated sum. units are dimensionless. - */ - doublereal sum_xlogx() const; - - //! Evaluate \f$ \sum_k X_k \log Q_k \f$. - /*! - * @param Q Vector of length m_kk to take the log average of - * @return Returns the indicated sum. - */ - doublereal sum_xlogQ(doublereal* const Q) const; - //@} - - /// @name Thermodynamic Properties - /// Class State only stores enough thermodynamic data to - /// specify the state. In addition to composition information, - /// it stores the temperature and - /// mass density. - //@{ - - //! Temperature (K). - /*! - * @return Returns the temperature of the phase - */ - doublereal temperature() const { - return m_temp; - } - - //! Density (kg/m^3). - /*! - * @return Returns the density of the phase - */ - virtual doublereal density() const { - return m_dens; - } - - //! Molar density (kmol/m^3). - /*! - * @return Returns the molar density of the phase - */ - doublereal molarDensity() const; - - //! Molar volume (m^3/kmol). - /*! - * @return Returns the molar volume of the phase - */ - doublereal molarVolume() const; - - //! Set the internally stored density (kg/m^3) of the phase - /*! - * Note the density of a phase is an indepedent variable. - * - * @param density Input density (kg/m^3). - */ - virtual void setDensity(const doublereal density) { - m_dens = density; - } - - //! Set the internally stored molar density (kmol/m^3) of the phase. - /*! - * @param molarDensity Input molar density (kmol/m^3). - */ - virtual void setMolarDensity(const doublereal molarDensity); - - //! Set the temperature (K). - /*! - * This function sets the internally stored temperature of the phase. - * - * @param temp Temperature in kelvin - */ - virtual void setTemperature(const doublereal temp) { - m_temp = temp; - } - //@} - - //! True if the number species has been set - bool ready() const; - - //! Every time the mole fractions have changed, this routine - //! will increment the stateMFNumber - /*! - * @param forceChange If this is true then the stateMFNumber always - * changes. This defaults to false. - */ - void stateMFChangeCalc(bool forceChange = false); - - //! Return the state number - int stateMFNumber() const; - -protected: - - /** - * @internal - * Initialize. Make a local copy of the vector of - * molecular weights, and resize the composition arrays to - * the appropriate size. The only information an instance of - * State has about the species is their molecular weights. - * - * @param mw Vector of molecular weights of the species. - */ - void init(const vector_fp& mw); //, density_is_independent = true); - - /** - * m_kk is the number of species in the phase - */ - size_t m_kk; - - //! Set the molecular weight of a single species to a given value - /*! - * @param k id of the species - * @param mw Molecular Weight (kg kmol-1) - */ - void setMolecularWeight(const int k, const double mw) { - m_molwts[k] = mw; - m_rmolwts[k] = 1.0/mw; - } - -private: - - /** - * Temperature. This is an independent variable - * units = Kelvin - */ - doublereal m_temp; - - /** - * Density. This is an independent variable except in - * the incompressible degenerate case. Thus, - * the pressure is determined from this variable - * not the other way round. - * units = kg m-3 - */ - doublereal m_dens; - - /** - * m_mmw is the mean molecular weight of the mixture - * (kg kmol-1) - */ - doublereal m_mmw; - - /** - * m_ym[k] = mole fraction of species k divided by the - * mean molecular weight of mixture. - */ - mutable vector_fp m_ym; - - /** - * m_y[k] = mass fraction of species k - */ - mutable vector_fp m_y; - - /** - * m_molwts[k] = molecular weight of species k (kg kmol-1) - */ - vector_fp m_molwts; - - /** - * m_rmolwts[k] = inverse of the molecular weight of species k - * units = kmol kg-1. - */ - vector_fp m_rmolwts; - - //! State Change variable - /*! - * Whenever the mole fraction vector changes, this int is - * incremented. - */ - int m_stateNum; - -}; - -//! Return the State Mole Fraction Number -inline int State::stateMFNumber() const -{ - return m_stateNum; -} - -} - -#endif diff --git a/include/cantera/thermo/SurfPhase.h b/include/cantera/thermo/SurfPhase.h index 6d3499e9d..6aca0f6ce 100644 --- a/include/cantera/thermo/SurfPhase.h +++ b/include/cantera/thermo/SurfPhase.h @@ -590,7 +590,7 @@ public: * in kmol/m2, using m_n0, the surface site density, * and size(k), which is defined to be the number of * surface sites occupied by the kth molecule. - * It then calls State::setConcentrations to set the + * It then calls Phase::setConcentrations to set the * internal concentration in the object. * * @param theta This is the surface site fraction @@ -607,7 +607,7 @@ public: * in kmol/m2, using m_n0, the surface site density, * and size(k), which is defined to be the number of * surface sites occupied by the kth molecule. - * It then calls State::setConcentrations to set the + * It then calls Phase::setConcentrations to set the * internal concentration in the object. * * @param theta This is the surface site fraction diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index 12b50d16c..d2025f6bc 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -327,7 +327,7 @@ void DebyeHuckel::setPressure(doublereal p) void DebyeHuckel::setState_TP(doublereal t, doublereal p) { - State::setTemperature(t); + Phase::setTemperature(t); /* * Store the current pressure */ @@ -387,7 +387,7 @@ void DebyeHuckel::calcDensity() vtotal += vbar[i] * x[i]; } doublereal dd = meanMolecularWeight() / vtotal; - State::setDensity(dd); + Phase::setDensity(dd); } diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index 5ffd87579..d5a043796 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -115,32 +115,32 @@ GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const void GibbsExcessVPSSTP::setMassFractions(const doublereal* const y) { - State::setMassFractions(y); + Phase::setMassFractions(y); getMoleFractions(DATA_PTR(moleFractions_)); } void GibbsExcessVPSSTP::setMassFractions_NoNorm(const doublereal* const y) { - State::setMassFractions_NoNorm(y); + Phase::setMassFractions_NoNorm(y); getMoleFractions(DATA_PTR(moleFractions_)); } void GibbsExcessVPSSTP::setMoleFractions(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); getMoleFractions(DATA_PTR(moleFractions_)); } void GibbsExcessVPSSTP::setMoleFractions_NoNorm(const doublereal* const x) { - State::setMoleFractions_NoNorm(x); + Phase::setMoleFractions_NoNorm(x); getMoleFractions(DATA_PTR(moleFractions_)); } void GibbsExcessVPSSTP::setConcentrations(const doublereal* const c) { - State::setConcentrations(c); + Phase::setConcentrations(c); getMoleFractions(DATA_PTR(moleFractions_)); } @@ -192,13 +192,13 @@ void GibbsExcessVPSSTP::calcDensity() vtotal += vbar[i] * moleFractions_[i]; } doublereal dd = meanMolecularWeight() / vtotal; - State::setDensity(dd); + Phase::setDensity(dd); delete [] vbar; } void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) { - State::setTemperature(t); + Phase::setTemperature(t); /* * Store the current pressure */ diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 5b11b86e4..14abde555 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -787,7 +787,7 @@ void HMWSoln::calcDensity() vtotal += vbar[i] * x[i]; } doublereal dd = meanMolecularWeight() / vtotal; - State::setDensity(dd); + Phase::setDensity(dd); } /* @@ -828,7 +828,7 @@ doublereal HMWSoln::thermalExpansionCoeff() const double HMWSoln::density() const { // calcDensity(); - return State::density(); + return Phase::density(); } /* @@ -893,7 +893,7 @@ void HMWSoln::setTemperature(const doublereal temp) */ void HMWSoln::setState_TP(doublereal temp, doublereal pres) { - State::setTemperature(temp); + Phase::setTemperature(temp); /* * Store the current pressure */ diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp index 45cafaeec..d58306511 100644 --- a/src/thermo/IdealMolalSoln.cpp +++ b/src/thermo/IdealMolalSoln.cpp @@ -304,7 +304,7 @@ void IdealMolalSoln::calcDensity() vtotal += vbar[i] * x[i]; } doublereal dd = meanMolecularWeight() / vtotal; - State::setDensity(dd); + Phase::setDensity(dd); } /* @@ -374,7 +374,7 @@ void IdealMolalSoln::setDensity(const doublereal rho) */ void IdealMolalSoln::setMolarDensity(const doublereal conc) { - double concI = State::molarDensity(); + double concI = Phase::molarDensity(); if (conc != concI) { throw CanteraError("IdealMolalSoln::setMolarDensity", "molarDensity/denisty is not an independent variable"); @@ -383,7 +383,7 @@ void IdealMolalSoln::setMolarDensity(const doublereal conc) void IdealMolalSoln::setState_TP(doublereal temp, doublereal pres) { - State::setTemperature(temp); + Phase::setTemperature(temp); m_Pcurrent = pres; updateStandardStateThermo(); //m_densWaterSS = m_waterSS->density(); diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp index 0227644c6..5d2b95f54 100644 --- a/src/thermo/IdealSolidSolnPhase.cpp +++ b/src/thermo/IdealSolidSolnPhase.cpp @@ -283,10 +283,10 @@ void IdealSolidSolnPhase::calcDensity() m_speciesMolarVolume.end(), dtmp); /* * Set the density in the parent State object directly, - * by calling the State::setDensity() function. + * by calling the Phase::setDensity() function. */ double dens = 1.0/invDens; - State::setDensity(dens); + Phase::setDensity(dens); } /** @@ -363,7 +363,7 @@ void IdealSolidSolnPhase::setMolarDensity(const doublereal n) */ void IdealSolidSolnPhase::setMoleFractions(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); calcDensity(); } @@ -374,7 +374,7 @@ void IdealSolidSolnPhase::setMoleFractions(const doublereal* const x) */ void IdealSolidSolnPhase::setMoleFractions_NoNorm(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); calcDensity(); } @@ -385,7 +385,7 @@ void IdealSolidSolnPhase::setMoleFractions_NoNorm(const doublereal* const x) */ void IdealSolidSolnPhase::setMassFractions(const doublereal* const y) { - State::setMassFractions(y); + Phase::setMassFractions(y); calcDensity(); } @@ -396,7 +396,7 @@ void IdealSolidSolnPhase::setMassFractions(const doublereal* const y) */ void IdealSolidSolnPhase::setMassFractions_NoNorm(const doublereal* const y) { - State::setMassFractions_NoNorm(y); + Phase::setMassFractions_NoNorm(y); calcDensity(); } @@ -407,7 +407,7 @@ void IdealSolidSolnPhase::setMassFractions_NoNorm(const doublereal* const y) */ void IdealSolidSolnPhase::setConcentrations(const doublereal* const c) { - State::setConcentrations(c); + Phase::setConcentrations(c); calcDensity(); } diff --git a/src/thermo/IdealSolnGasVPSS.cpp b/src/thermo/IdealSolnGasVPSS.cpp index f713b19de..ae7092113 100644 --- a/src/thermo/IdealSolnGasVPSS.cpp +++ b/src/thermo/IdealSolnGasVPSS.cpp @@ -185,17 +185,17 @@ void IdealSolnGasVPSS::calcDensity() if (m_idealGas) { double dens = (m_Pcurrent * meanMolecularWeight() /(GasConstant * temperature())); - State::setDensity(dens); + Phase::setDensity(dens); } else { const doublereal* const dtmp = moleFractdivMMW(); const vector_fp& vss = m_VPSS_ptr->standardVolumes(); double invDens = dot(vss.begin(), vss.end(), dtmp); /* * Set the density in the parent State object directly, - * by calling the State::setDensity() function. + * by calling the Phase::setDensity() function. */ double dens = 1.0/invDens; - State::setDensity(dens); + Phase::setDensity(dens); } } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index fd05bfc5f..10d7039aa 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -629,7 +629,7 @@ void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p) //calcDensity(); double dd = neutralMoleculePhase_->density(); - State::setDensity(dd); + Phase::setDensity(dd); } // Calculate ion mole fractions from neutral molecule diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index d86b774c4..fc268fdc5 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -239,10 +239,10 @@ doublereal LatticePhase::calcDensity() // doublereal invDens = dot(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), dtmp); /* * Set the density in the parent State object directly, - * by calling the State::setDensity() function. + * by calling the Phase::setDensity() function. */ // doublereal dens = 1.0/invDens; - // State::setDensity(dens); + // Phase::setDensity(dens); return dens; } //==================================================================================================================== @@ -254,31 +254,31 @@ void LatticePhase::setPressure(doublereal p) //==================================================================================================================== void LatticePhase::setMoleFractions(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); calcDensity(); } //==================================================================================================================== void LatticePhase::setMoleFractions_NoNorm(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); calcDensity(); } //==================================================================================================================== void LatticePhase::setMassFractions(const doublereal* const y) { - State::setMassFractions(y); + Phase::setMassFractions(y); calcDensity(); } //==================================================================================================================== void LatticePhase::setMassFractions_NoNorm(const doublereal* const y) { - State::setMassFractions_NoNorm(y); + Phase::setMassFractions_NoNorm(y); calcDensity(); } //==================================================================================================================== void LatticePhase::setConcentrations(const doublereal* const c) { - State::setConcentrations(c); + Phase::setConcentrations(c); calcDensity(); } //==================================================================================================================== diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index 2a3bc2cda..fd58500fc 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -284,7 +284,7 @@ doublereal LatticeSolidPhase::calcDensity() for (size_t n = 0; n < m_nlattice; n++) { sum += theta_[n] * m_lattice[n]->density(); } - State::setDensity(sum); + Phase::setDensity(sum); return sum; } //==================================================================================================================== @@ -312,7 +312,7 @@ void LatticeSolidPhase::setMoleFractions(const doublereal* const x) for (size_t k = 0; k < strt; k++) { m_x[k] = x[k] / m_nlattice; } - State::setMoleFractions(DATA_PTR(m_x)); + Phase::setMoleFractions(DATA_PTR(m_x)); calcDensity(); } //==================================================================================================================== @@ -327,7 +327,7 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const { size_t nsp, strt = 0; // the ifdef block should be the way we calculate this.!!!!! - State::getMoleFractions(x); + Phase::getMoleFractions(x); doublereal sum; for (size_t n = 0; n < m_nlattice; n++) { nsp = m_lattice[n]->nSpecies(); diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index 15fb0ab0d..b2187c9c0 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -476,7 +476,7 @@ void MixtureFugacityTP::setStateFromXML(const XML_Node& state) double rho = ctml::getFloat(state, "density", "density"); setState_TR(t, rho); } else if (doTP) { - double rho = State::density(); + double rho = Phase::density(); setState_TR(t, rho); } } @@ -525,37 +525,37 @@ void MixtureFugacityTP::setPressure(doublereal p) //==================================================================================================================== void MixtureFugacityTP::setMassFractions(const doublereal* const y) { - State::setMassFractions(y); + Phase::setMassFractions(y); getMoleFractions(DATA_PTR(moleFractions_)); } //==================================================================================================================== void MixtureFugacityTP::setMassFractions_NoNorm(const doublereal* const y) { - State::setMassFractions_NoNorm(y); + Phase::setMassFractions_NoNorm(y); getMoleFractions(DATA_PTR(moleFractions_)); } //==================================================================================================================== void MixtureFugacityTP::setMoleFractions(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); getMoleFractions(DATA_PTR(moleFractions_)); } //==================================================================================================================== void MixtureFugacityTP::setMoleFractions_NoNorm(const doublereal* const x) { - State::setMoleFractions_NoNorm(x); + Phase::setMoleFractions_NoNorm(x); getMoleFractions(DATA_PTR(moleFractions_)); } //==================================================================================================================== void MixtureFugacityTP::setConcentrations(const doublereal* const c) { - State::setConcentrations(c); + Phase::setConcentrations(c); getMoleFractions(DATA_PTR(moleFractions_)); } //==================================================================================================================== void MixtureFugacityTP::setMoleFractions_NoState(const doublereal* const x) { - State::setMoleFractions(x); + Phase::setMoleFractions(x); getMoleFractions(DATA_PTR(moleFractions_)); updateMixingExpressions(); } @@ -580,7 +580,7 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) getMoleFractions(DATA_PTR(moleFractions_)); - State::setTemperature(t); + Phase::setTemperature(t); _updateReferenceStateThermo(); // Depends on the mole fractions and the temperature updateMixingExpressions(); @@ -589,17 +589,17 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) // double mmw = meanMolecularWeight(); if (forcedState_ == FLUID_UNDEFINED) { - double rhoNow = State::density(); + double rhoNow = Phase::density(); double rho = densityCalc(t, pres, iState_, rhoNow); if (rho > 0.0) { - State::setDensity(rho); + Phase::setDensity(rho); m_Pcurrent = pres; iState_ = phaseState(true); } else { if (rho < -1.5) { rho = densityCalc(t, pres, FLUID_UNDEFINED , rhoNow); if (rho > 0.0) { - State::setDensity(rho); + Phase::setDensity(rho); m_Pcurrent = pres; iState_ = phaseState(true); } else { @@ -615,10 +615,10 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) } else if (forcedState_ == FLUID_GAS) { // Normal density calculation if (iState_ < FLUID_LIQUID_0) { - double rhoNow = State::density(); + double rhoNow = Phase::density(); double rho = densityCalc(t, pres, iState_, rhoNow); if (rho > 0.0) { - State::setDensity(rho); + Phase::setDensity(rho); m_Pcurrent = pres; iState_ = phaseState(true); if (iState_ >= FLUID_LIQUID_0) { @@ -633,10 +633,10 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) } else if (forcedState_ > FLUID_LIQUID_0) { if (iState_ >= FLUID_LIQUID_0) { - double rhoNow = State::density(); + double rhoNow = Phase::density(); double rho = densityCalc(t, pres, iState_, rhoNow); if (rho > 0.0) { - State::setDensity(rho); + Phase::setDensity(rho); m_Pcurrent = pres; iState_ = phaseState(true); if (iState_ == FLUID_GAS) { @@ -667,9 +667,9 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho) { getMoleFractions(DATA_PTR(moleFractions_)); - State::setTemperature(T); + Phase::setTemperature(T); _updateReferenceStateThermo(); - State::setDensity(rho); + Phase::setDensity(rho); doublereal mv = molarVolume(); // depends on mole fraction and temperature updateMixingExpressions(); diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp index 4956ddbbe..0a7542601 100644 --- a/src/thermo/Phase.cpp +++ b/src/thermo/Phase.cpp @@ -19,12 +19,15 @@ namespace Cantera Phase::Phase() : Constituents(), - State(), m_kk(0), m_ndim(3), m_xml(new XML_Node("phase")), m_id(""), - m_name("") + m_name(""), + m_temp(0.0), + m_dens(0.001), + m_mmw(0.0), + m_stateNum(-1) { } @@ -36,12 +39,15 @@ Phase::Phase() : */ Phase::Phase(const Phase& right) : Constituents(), - State(), m_kk(0), m_ndim(3), m_xml(new XML_Node("phase")), m_id(""), - m_name("") + m_name(""), + m_temp(0.0), + m_dens(0.001), + m_mmw(0.0), + m_stateNum(-1) { /* * Call the assignment operator. @@ -70,13 +76,21 @@ Phase& Phase::operator=(const Phase& right) * Now call the inherited-classes assignment operators. */ (void) Constituents::operator=(right); - (void) State::operator=(right); /* * Handle its own data */ m_kk = right.m_kk; m_ndim = right.m_ndim; m_data = right.m_data; + m_temp = right.m_temp; + m_dens = right.m_dens; + m_mmw = right.m_mmw; + m_ym = right.m_ym; + m_y = right.m_y; + m_molwts = right.m_molwts; + m_rmolwts = right.m_rmolwts; + m_stateNum = -1; + /* * This is a little complicated. -> Because we delete m_xml * in the destructor, we own m_xml completely, and we need @@ -106,6 +120,16 @@ Phase::~Phase() } } +inline void Phase::stateMFChangeCalc(bool forcerChange) +{ + // Right now we assume that the mole fractions have changed every time + // the function is called + m_stateNum++; + if (m_stateNum > 1000000) { + m_stateNum = -10000000; + } +} + XML_Node& Phase::xml() { return *m_xml; @@ -197,6 +221,32 @@ void Phase::restoreState(size_t lenstate, const doublereal* state) } } +void Phase::setMoleFractions(const doublereal* const x) +{ + doublereal sum = dot(x, x + m_kk, m_molwts.begin()); + doublereal rsum = 1.0/sum; + transform(x, x + m_kk, m_ym.begin(), timesConstant(rsum)); + transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), + m_y.begin(), multiplies()); + doublereal norm = accumulate(x, x + m_kk, 0.0); + m_mmw = sum/norm; + + // Call a routine to determine whether state has changed. + stateMFChangeCalc(); +} + +void Phase::setMoleFractions_NoNorm(const doublereal* const x) +{ + m_mmw = dot(x, x + m_kk, m_molwts.begin()); + doublereal rmmw = 1.0/m_mmw; + transform(x, x + m_kk, m_ym.begin(), timesConstant(rmmw)); + transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), + m_y.begin(), multiplies()); + + // Call a routine to determine whether state has changed. + stateMFChangeCalc(); +} + void Phase::setMoleFractionsByName(compositionMap& xMap) { size_t kk = nSpecies(); @@ -220,12 +270,35 @@ void Phase::setMoleFractionsByName(const std::string& x) } parseCompString(x, xx); setMoleFractionsByName(xx); - //int kk = nSpecies(); - //vector_fp mf(kk); - //for (int k = 0; k < kk; k++) { - // mf[k] = xx[speciesName(k)]; - //} - //setMoleFractions(mf.begin()); +} + +void Phase::setMassFractions(const doublereal* const y) +{ + doublereal norm = 0.0, sum = 0.0; + norm = accumulate(y, y + m_kk, 0.0); + copy(y, y + m_kk, m_y.begin()); + scale(y, y + m_kk, m_y.begin(), 1.0/norm); + + transform(m_y.begin(), m_y.begin() + m_kk, m_rmolwts.begin(), + m_ym.begin(), multiplies()); + sum = accumulate(m_ym.begin(), m_ym.begin() + m_kk, 0.0); + m_mmw = 1.0/sum; + + // Call a routine to determine whether state has changed. + stateMFChangeCalc(); +} + +void Phase::setMassFractions_NoNorm(const doublereal* const y) +{ + doublereal sum = 0.0; + copy(y, y + m_kk, m_y.begin()); + transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(), + multiplies()); + sum = accumulate(m_ym.begin(), m_ym.end(), 0.0); + m_mmw = 1.0/sum; + + // Call a routine to determine whether state has changed. + stateMFChangeCalc(); } void Phase::setMassFractionsByName(compositionMap& yMap) @@ -381,28 +454,48 @@ void Phase::getMoleFractionsByName(compositionMap& x) const x.clear(); size_t kk = nSpecies(); for (size_t k = 0; k < kk; k++) { - x[speciesName(k)] = State::moleFraction(k); + x[speciesName(k)] = Phase::moleFraction(k); } } +void Phase::getMoleFractions(doublereal* const x) const +{ + scale(m_ym.begin(), m_ym.end(), x, m_mmw); +} + doublereal Phase::moleFraction(size_t k) const { - return State::moleFraction(k); + if (k < m_kk) { + return m_ym[k] * m_mmw; + } else { + throw CanteraError("Phase::moleFraction", + "illegal species index number"); + } + return 0.0; } doublereal Phase::moleFraction(std::string nameSpec) const { size_t iloc = speciesIndex(nameSpec); if (iloc != npos) { - return State::moleFraction(iloc); + return moleFraction(iloc); } else { return 0.0; } } +const doublereal* Phase::moleFractdivMMW() const +{ + return &m_ym[0]; +} + doublereal Phase::massFraction(size_t k) const { - return State::massFraction(k); + if (k < m_kk) { + return m_y[k]; + } + throw CanteraError("State:massFraction", "illegal species index number"); + return 0.0; } doublereal Phase::massFraction(std::string nameSpec) const @@ -415,17 +508,90 @@ doublereal Phase::massFraction(std::string nameSpec) const } } +void Phase::getMassFractions(doublereal* const y) const +{ + copy(m_y.begin(), m_y.end(), y); +} + +doublereal Phase::concentration(const size_t k) const +{ + if (k < m_kk) { + return m_y[k] * m_dens * m_rmolwts[k] ; + } + throw CanteraError("State:massFraction", "illegal species index number"); + return 0.0; +} + +void Phase::getConcentrations(doublereal* const c) const +{ + scale(m_ym.begin(), m_ym.end(), c, m_dens); +} + +void Phase::setConcentrations(const doublereal* const conc) +{ + doublereal sum = 0.0, norm = 0.0; + for (size_t k = 0; k != m_kk; ++k) { + sum += conc[k]*m_molwts[k]; + norm += conc[k]; + } + m_mmw = sum/norm; + setDensity(sum); + doublereal rsum = 1.0/sum; + for (size_t k = 0; k != m_kk; ++k) { + m_ym[k] = conc[k] * rsum; + m_y[k] = m_ym[k] * m_molwts[k]; + } + + // Call a routine to determine whether state has changed. + stateMFChangeCalc(); +} + +doublereal Phase::molarDensity() const +{ + return density()/meanMolecularWeight(); +} + +void Phase::setMolarDensity(const doublereal molarDensity) +{ + m_dens = molarDensity*meanMolecularWeight(); +} + +doublereal Phase::molarVolume() const +{ + return 1.0/molarDensity(); +} + doublereal Phase::chargeDensity() const { size_t kk = nSpecies(); doublereal cdens = 0.0; for (size_t k = 0; k < kk; k++) { - cdens += charge(k)*State::moleFraction(k); + cdens += charge(k)*moleFraction(k); } cdens *= Faraday; return cdens; } +doublereal Phase::mean_X(const doublereal* const Q) const +{ + return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0); +} + +doublereal Phase::mean_Y(const doublereal* const Q) const +{ + return dot(m_y.begin(), m_y.end(), Q); +} + +doublereal Phase::sum_xlogx() const +{ + return m_mmw* Cantera::sum_xlogx(m_ym.begin(), m_ym.end()) + log(m_mmw); +} + +doublereal Phase::sum_xlogQ(doublereal* Q) const +{ + return m_mmw * Cantera::sum_xlogQ(m_ym.begin(), m_ym.end(), Q); +} + /** * Finished adding species, prepare to use them for calculation * of mixture properties. @@ -444,16 +610,47 @@ void Phase::freezeSpecies() m_kk = nSpecies(); } +void Phase::init(const vector_fp& mw) +{ + m_kk = mw.size(); + m_molwts.resize(m_kk); + m_rmolwts.resize(m_kk); + m_y.resize(m_kk, 0.0); + m_ym.resize(m_kk, 0.0); + copy(mw.begin(), mw.end(), m_molwts.begin()); + for (size_t k = 0; k < m_kk; k++) { + if (m_molwts[k] < 0.0) { + throw CanteraError("Phase::init", + "negative molecular weight for species number " + + int2str(k)); + } + /* + * Some surface phases may define species representing + * empty sites that have zero molecular weight. Give them + * a very small molecular weight to avoid dividing by + * zero. + */ + if (m_molwts[k] < Tiny) { + m_molwts[k] = Tiny; + } + m_rmolwts[k] = 1.0/m_molwts[k]; + } + + /* + * Now that we have resized the State object, let's fill it with + * a valid mass fraction vector that sums to one. The State object + * should never have a mass fraction vector that doesn't sum to one. + * We will assume that species 0 has a mass fraction of 1.0 and + * mass fraction of all other species is 0.0. + */ + m_y[0] = 1.0; + m_ym[0] = m_y[0] * m_rmolwts[0]; + m_mmw = 1.0 / m_ym[0]; +} + bool Phase::ready() const { - return (m_kk > 0 && Constituents::ready() && State::ready()); + return (m_kk > 0 && Constituents::ready()); } -// int Phase::installUpdater_T(Updater* u) { -// return m_T_updater.install(u); -// } - -// int Phase::installUpdater_C(Updater* u) { -// return m_C_updater.install(u); -// } -} +} // namespace Cantera diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 2271f4050..6b4c3b30c 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -371,17 +371,17 @@ void RedlichKwongMFTP::calcDensity() double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp); /* * Set the density in the parent State object directly, - * by calling the State::setDensity() function. + * by calling the Phase::setDensity() function. */ double dens = 1.0/invDens; - State::setDensity(dens); + Phase::setDensity(dens); } //==================================================================================================================== void RedlichKwongMFTP::setTemperature(const doublereal temp) { - State::setTemperature(temp); + Phase::setTemperature(temp); _updateReferenceStateThermo(); updateAB(); } diff --git a/src/thermo/State.cpp b/src/thermo/State.cpp deleted file mode 100644 index fa0fd98e8..000000000 --- a/src/thermo/State.cpp +++ /dev/null @@ -1,311 +0,0 @@ -/** - * @file State.cpp - * Definitions for the class State, that manages the independent variables of temperature, mass density, - * and species mass/mole fraction that define the thermodynamic state (see \ref phases and - * class \link Cantera::State State\endlink). - */ - -/* - * Copyright 2003-2004 California Institute of Technology - * See file License.txt for licensing information - */ - -#include "cantera/base/utilities.h" -#include "cantera/base/ctexceptions.h" -#include "cantera/base/stringUtils.h" -#include "cantera/thermo/State.h" - -//#ifdef DARWIN -//#include -//#endif - -using namespace std; - -namespace Cantera -{ - -inline void State::stateMFChangeCalc(bool forcerChange) -{ - // Right now we assume that the mole fractions have changed every time - // the function is called - m_stateNum++; - if (m_stateNum > 1000000) { - m_stateNum = -10000000; - } -} - -State::State() : - m_kk(0), - m_temp(0.0), - m_dens(0.001), - m_mmw(0.0), - m_stateNum(-1) -{ -} - -State::~State() -{ -} - -State::State(const State& right) : - m_kk(0), - m_temp(0.0), - m_dens(0.001), - m_mmw(0.0), - m_stateNum(-1) -{ - /* - * Call the assignment operator. - */ - *this = operator=(right); -} - -/* - * Assignment operator for the State Class - */ -State& State::operator=(const State& right) -{ - /* - * Check for self assignment. - */ - if (this == &right) { - return *this; - } - /* - * We do a straight assignment operator on all of the - * data. The vectors are copied. - */ - m_kk = right.m_kk; - m_temp = right.m_temp; - m_dens = right.m_dens; - m_mmw = right.m_mmw; - m_ym = right.m_ym; - m_y = right.m_y; - m_molwts = right.m_molwts; - m_rmolwts = right.m_rmolwts; - m_stateNum = -1; - /* - * Return the reference to the current object - */ - return *this; -} - -doublereal State::moleFraction(const size_t k) const -{ - if (k < m_kk) { - return m_ym[k] * m_mmw; - } else { - throw CanteraError("State:moleFraction", - "illegal species index number"); - } - return 0.0; -} - -void State::setMoleFractions(const doublereal* const x) -{ - doublereal sum = dot(x, x + m_kk, m_molwts.begin()); - doublereal rsum = 1.0/sum; - transform(x, x + m_kk, m_ym.begin(), timesConstant(rsum)); - transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), - m_y.begin(), multiplies()); - doublereal norm = accumulate(x, x + m_kk, 0.0); - m_mmw = sum/norm; - - //! Call a routine to determine whether state has changed. - stateMFChangeCalc(); -} - -void State::setMoleFractions_NoNorm(const doublereal* const x) -{ - m_mmw = dot(x, x + m_kk, m_molwts.begin()); - doublereal rmmw = 1.0/m_mmw; - transform(x, x + m_kk, m_ym.begin(), timesConstant(rmmw)); - transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), - m_y.begin(), multiplies()); - - //! Call a routine to determine whether state has changed. - stateMFChangeCalc(); -} - -doublereal State::massFraction(const size_t k) const -{ - if (k < m_kk) { - return m_y[k]; - } - throw CanteraError("State:massFraction", "illegal species index number"); - return 0.0; -} - -doublereal State::concentration(const size_t k) const -{ - if (k < m_kk) { - return m_y[k] * m_dens * m_rmolwts[k] ; - } - throw CanteraError("State:massFraction", "illegal species index number"); - return 0.0; -} - -void State::setMassFractions(const doublereal* const y) -{ - doublereal norm = 0.0, sum = 0.0; - //cblas_dcopy(m_kk, y, 1, m_y.begin(), 1); - norm = accumulate(y, y + m_kk, 0.0); - copy(y, y + m_kk, m_y.begin()); - scale(y, y + m_kk, m_y.begin(), 1.0/norm); - // for (k = 0; k != m_kk; ++k) { - // norm += y[k]; - // m_y[k] = y[k]; - //} - - //scale(m_kk, 1.0/norm, m_y.begin()); - transform(m_y.begin(), m_y.begin() + m_kk, m_rmolwts.begin(), - m_ym.begin(), multiplies()); - sum = accumulate(m_ym.begin(), m_ym.begin() + m_kk, 0.0); - // for (k = 0; k != m_kk; ++k) { - // m_ym[k] = m_y[k] * m_rmolwts[k]; - // sum += m_ym[k]; - // } - m_mmw = 1.0/sum; - - //! Call a routine to determine whether state has changed. - stateMFChangeCalc(); -} - -void State::setMassFractions_NoNorm(const doublereal* const y) -{ - doublereal sum = 0.0; - copy(y, y + m_kk, m_y.begin()); - transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(), - multiplies()); - sum = accumulate(m_ym.begin(), m_ym.end(), 0.0); - //for (k = 0; k != m_kk; ++k) { - // m_y[k] = y[k]; - // m_ym[k] = m_y[k] * m_rmolwts[k]; - // sum += m_ym[k]; - //} - m_mmw = 1.0/sum; - - //! Call a routine to determine whether state has changed. - stateMFChangeCalc(); -} - -doublereal State::sum_xlogx() const -{ - return m_mmw* Cantera::sum_xlogx(m_ym.begin(), m_ym.end()) + log(m_mmw); -} - -doublereal State::sum_xlogQ(doublereal* Q) const -{ - return m_mmw * Cantera::sum_xlogQ(m_ym.begin(), m_ym.end(), Q); -} - -doublereal State::molarDensity() const -{ - return density()/meanMolecularWeight(); -} - -doublereal State::molarVolume() const -{ - return 1.0/molarDensity(); -} - -void State::setConcentrations(const doublereal* const conc) -{ - doublereal sum = 0.0, norm = 0.0; - for (size_t k = 0; k != m_kk; ++k) { - sum += conc[k]*m_molwts[k]; - norm += conc[k]; - } - m_mmw = sum/norm; - setDensity(sum); - doublereal rsum = 1.0/sum; - for (size_t k = 0; k != m_kk; ++k) { - m_ym[k] = conc[k] * rsum; - m_y[k] = m_ym[k] * m_molwts[k]; - } - - // Call a routine to determine whether state has changed. - stateMFChangeCalc(); -} - -const doublereal* State::moleFractdivMMW() const -{ - return &m_ym[0]; -} - -void State::getConcentrations(doublereal* const c) const -{ - scale(m_ym.begin(), m_ym.end(), c, m_dens); -} - -doublereal State::mean_X(const doublereal* const Q) const -{ - return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0); -} - -doublereal State::mean_Y(const doublereal* const Q) const -{ - return dot(m_y.begin(), m_y.end(), Q); -} - -void State::getMoleFractions(doublereal* const x) const -{ - scale(m_ym.begin(), m_ym.end(), x, m_mmw); -} - -void State::getMassFractions(doublereal* const y) const -{ - copy(m_y.begin(), m_y.end(), y); -} - -void State::setMolarDensity(const doublereal molarDensity) -{ - m_dens = molarDensity*meanMolecularWeight(); -} - -void State::init(const vector_fp& mw) -{ - m_kk = mw.size(); - m_molwts.resize(m_kk); - m_rmolwts.resize(m_kk); - m_y.resize(m_kk, 0.0); - m_ym.resize(m_kk, 0.0); - copy(mw.begin(), mw.end(), m_molwts.begin()); - for (size_t k = 0; k < m_kk; k++) { - if (m_molwts[k] < 0.0) { - throw CanteraError("State::init", - "negative molecular weight for species number " - + int2str(k)); - } - /* - * Some surface phases may define species representing - * empty sites that have zero molecular weight. Give them - * a very small molecular weight to avoid dividing by - * zero. - */ - if (m_molwts[k] < Tiny) { - m_molwts[k] = Tiny; - } - m_rmolwts[k] = 1.0/m_molwts[k]; - } - - /* - * Now that we have resized the State object, let's fill it with - * a valid mass fraction vector that sums to one. The State object - * should never have a mass fraction vector that doesn't sum to one. - * We will assume that species 0 has a mass fraction of 1.0 and - * mass fraction of all other species is 0.0. - */ - m_y[0] = 1.0; - m_ym[0] = m_y[0] * m_rmolwts[0]; - m_mmw = 1.0 / m_ym[0]; -} - -// True if the number of species has been set and fixed -bool State::ready() const -{ - return (m_kk > 0); -} - - -} diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp index cc68720c7..8dc5864c1 100644 --- a/src/thermo/SurfPhase.cpp +++ b/src/thermo/SurfPhase.cpp @@ -375,7 +375,7 @@ void SurfPhase::setSiteDensity(doublereal n0) * in kmol/m2, using m_n0, the surface site density, * and size(k), which is defined to be the number of * surface sites occupied by the kth molecule. - * It then calls State::setConcentrations to set the + * It then calls Phase::setConcentrations to set the * internal concentration in the object. */ void SurfPhase:: @@ -396,7 +396,7 @@ setCoverages(const doublereal* theta) m_work[k] = m_n0*theta[k]/(sum*size(k)); } /* - * Call the State:: class function + * Call the Phase:: class function * setConcentrations. */ setConcentrations(DATA_PTR(m_work)); @@ -409,7 +409,7 @@ setCoveragesNoNorm(const doublereal* theta) m_work[k] = m_n0*theta[k]/(size(k)); } /* - * Call the State:: class function + * Call the Phase:: class function * setConcentrations. */ setConcentrations(DATA_PTR(m_work)); diff --git a/src/thermo/VPStandardStateTP.cpp b/src/thermo/VPStandardStateTP.cpp index d97be06cf..aa4767ebc 100644 --- a/src/thermo/VPStandardStateTP.cpp +++ b/src/thermo/VPStandardStateTP.cpp @@ -417,7 +417,7 @@ void VPStandardStateTP::setState_TP(doublereal t, doublereal pres) * Therefore, we need to do the standard state thermo calc with the * (t, pres) combo. */ - State::setTemperature(t); + Phase::setTemperature(t); m_Pcurrent = pres; updateStandardStateThermo(); /* diff --git a/src/thermo/WaterSSTP.cpp b/src/thermo/WaterSSTP.cpp index ff48698f1..bd7566211 100644 --- a/src/thermo/WaterSSTP.cpp +++ b/src/thermo/WaterSSTP.cpp @@ -246,8 +246,8 @@ initThermoXML(XML_Node& phaseNode, std::string id) * Set the baseline */ doublereal T = 298.15; - State::setDensity(7.0E-8); - State::setTemperature(T); + Phase::setDensity(7.0E-8); + Phase::setTemperature(T); doublereal presLow = 1.0E-2; doublereal oneBar = 1.0E5; @@ -588,14 +588,14 @@ doublereal WaterSSTP::critDensity() const void WaterSSTP::setTemperature(const doublereal temp) { - State::setTemperature(temp); + Phase::setTemperature(temp); doublereal dd = density(); m_sub->setState_TR(temp, dd); } void WaterSSTP::setDensity(const doublereal dens) { - State::setDensity(dens); + Phase::setDensity(dens); doublereal temp = temperature(); m_sub->setState_TR(temp, dens); }