From 1e07d47bb86830658492b277395b8d66b738fb5e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 28 Mar 2007 20:33:43 +0000 Subject: [PATCH] Added an element abundance return capability. doxygen update Started documenting MultiPhase --- Cantera/src/MultiPhase.cpp | 53 ++++ Cantera/src/MultiPhase.h | 585 ++++++++++++++++++++----------------- tools/doc/Cantera.cfg.in | 2 +- 3 files changed, 371 insertions(+), 269 deletions(-) diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index 36ed50690..305a39d66 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -118,6 +118,8 @@ namespace Cantera { m_moleFractions.resize(m_nsp, 0.0); // iterate over the elements + // -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), + // m_sptart(ip) for (m = 0; m < m_nel; m++) { sym = m_enames[m]; k = 0; @@ -158,6 +160,10 @@ namespace Cantera { /// mole fractions stored in the phase objects m_init = true; updateMoleFractions(); + + m_elemAbundances.resize(m_nel, 0.0); + + calcElemAbundances(); } @@ -415,6 +421,38 @@ namespace Cantera { } } + void MultiPhase::getElemAbundances(doublereal *elemAbundances) const { + index_t eGlobal; + calcElemAbundances(); + for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + elemAbundances[eGlobal] = m_elemAbundances[eGlobal]; + } + } + + // Internal routine to calculate the element abundance vector + void MultiPhase::calcElemAbundances() const { + index_t loc = 0; + index_t eGlobal; + int ik, kGlobal; + doublereal spMoles; + for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + m_elemAbundances[eGlobal] = 0.0; + } + for (index_t ip = 0; ip < m_np; ip++) { + phase_t* p = m_phase[ip]; + int nspPhase = p->nSpecies(); + doublereal phasemoles = m_moles[ip]; + loc += nspPhase; + for (ik = 0; ik < nspPhase; ik++) { + kGlobal = loc + ik; + spMoles = m_moleFractions[kGlobal] * phasemoles; + for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles; + } + } + } + } + /// The total mixture volume [m^3]. doublereal MultiPhase::volume() const { int i; @@ -772,6 +810,21 @@ done: } #endif + // Name of element \a m. + std::string MultiPhase::elementName(int m) const { + return m_enames[m]; + } + + // Index of element with name \a name. + int MultiPhase::elementIndex(std::string name) const { + return m_enamemap[name] - 1; + } + + // Name of species with global index \a k. + std::string MultiPhase::speciesName(int k) const { + return m_snames[k]; + } + //------------------------------------------------------------- // // protected methods diff --git a/Cantera/src/MultiPhase.h b/Cantera/src/MultiPhase.h index b81901157..94c3fb07e 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -1,6 +1,10 @@ /** * @file MultiPhase.h + * Headers for the \link Cantera::MultiPhase MultiPhase\endlink + * object that is used to set up multiphase equilibrium problems (see \ref equilfunctions). * + */ +/* * $Author$ * $Date$ * $Revision$ @@ -14,290 +18,335 @@ namespace Cantera { - /// A class for multiphase mixtures. The mixture can contain any - /// number of phases of any type. All phases have the same - /// temperature and pressure, and a specified number of moles. - /// The phases do not need to have the same elements. For example, - /// a mixture might consist of a gaseous phase with elements (H, - /// C, O, N), a solid carbon phase containing only element C, - /// etc. A master element set will be constructed for the mixture - /// that is the union of the elements of each phase. + //! A class for multiphase mixtures. The mixture can contain any + //! number of phases of any type. + /*! + * All phases have the same + * temperature and pressure, and a specified number of moles. + * The phases do not need to have the same elements. For example, + * a mixture might consist of a gaseous phase with elements (H, + * C, O, N), a solid carbon phase containing only element C, + * etc. A master element set will be constructed for the mixture + * that is the union of the elements of each phase. + * + * @ingroup equilfunctions + */ + class MultiPhase { - class MultiPhase { + public: - public: + // some typedefs for convenience + typedef size_t index_t; + typedef ThermoPhase phase_t; + typedef DenseMatrix array_t; + typedef std::vector phase_list; - // some typedefs for convenience - typedef size_t index_t; - typedef ThermoPhase phase_t; - typedef DenseMatrix array_t; - typedef std::vector phase_list; + /// Constructor. The constructor takes no arguments, since + /// phases are added using method addPhase. + MultiPhase(); - /// Constructor. The constructor takes no arguments, since - /// phases are added using method addPhase. - MultiPhase(); - - /// Destructor. Does nothing. Class MultiPhase does not take - /// "ownership" (i.e. responsibility for destroying) the - /// phase objects. - virtual ~MultiPhase() {} + /// Destructor. Does nothing. Class MultiPhase does not take + /// "ownership" (i.e. responsibility for destroying) the + /// phase objects. + virtual ~MultiPhase() {} - void addPhases(phase_list& phases, const vector_fp& phaseMoles); + void addPhases(phase_list& phases, const vector_fp& phaseMoles); - /// Add all phases present in 'mix' to this mixture. - void addPhases(MultiPhase& mix); + /// Add all phases present in 'mix' to this mixture. + void addPhases(MultiPhase& mix); - /// Add a phase to the mixture. - /// @param p pointer to the phase object - /// @param moles total number of moles of all species in this phase - void addPhase(phase_t* p, doublereal moles); + /// Add a phase to the mixture. + /// @param p pointer to the phase object + /// @param moles total number of moles of all species in this phase + void addPhase(phase_t* p, doublereal moles); - /// Number of elements. - int nElements() const { return int(m_nel); } + /// Number of elements. + int nElements() const { return int(m_nel); } - /// Name of element \a m. - std::string elementName(int m) const { return m_enames[m]; } + /// Name of element \a m. + std::string elementName(int m) const; - /// Index of element with name \a name. - int elementIndex(std::string name) const { return m_enamemap[name] - 1;} + /// Index of element with name \a name. + int elementIndex(std::string name) const; - /// Number of species, summed over all phases. - int nSpecies() const { return int(m_nsp); } + //! Number of species, summed over all phases. + int nSpecies() const { return int(m_nsp); } - /// Name of species with index \a k. - std::string speciesName(int k) const { return m_snames[k]; } + //! Name of species with global index \a k. + std::string speciesName(int k) const; - /// Number of atoms of element \a m in species \a k. - doublereal nAtoms(int k, int m) { - if (!m_init) init(); - return m_atoms(m,k); - } - - /// Species mole fractions. Write the array of species mole - /// fractions into array \c x. The mole fractions are - /// normalized to sum to one in each phase. - void getMoleFractions(doublereal* x) const { - std::copy(m_moleFractions.begin(), m_moleFractions.end(), x); - } - - /// Process phases and build atomic composition array. After - /// init() has been called, no more phases may be added. - void init(); - - /// Moles of phase n. - doublereal phaseMoles(index_t n) const { - return m_moles[n]; - } - - /// Set the number of moles of phase with index n. - void setPhaseMoles(index_t n, doublereal moles) { - m_moles[n] = moles; - } - - /// Return a reference to phase n. The state of phase n is - /// also updated to match the state stored locally in the - /// mixture object. - phase_t& phase(index_t n); - - /// Moles of species \c k. - doublereal speciesMoles(index_t k) const; - - /// Index of the species belonging to phase number \c p - /// with index \c k within the phase. - int speciesIndex(index_t k, index_t p) const { - return m_spstart[p] + k; - } - - /// Minimum temperature for which all solution phases have - /// valid thermo data. Stoichiometric phases are not - /// considered, since they may have thermo data only valid for - /// conditions for which they are stable. - doublereal minTemp() const { return m_Tmin; } - - /// Maximum temperature for which all solution phases have - /// valid thermo data. Stoichiometric phases are not - /// considered, since they may have thermo data only valid for - /// conditions for which they are stable. - doublereal maxTemp() const { return m_Tmax; } - - /// Total charge (Coulombs). - doublereal charge() const; - - /// Charge (Coulombs) of phase with index \a p. - doublereal phaseCharge(index_t p) const; - - /// Total moles of element \a m, summed over all phases. - doublereal elementMoles(index_t m) const; - - /// Chemical potentials. Write into array \a mu the chemical - /// potentials of all species [J/kmol]. The chemical - /// potentials are related to the activities by - /// \f[ \mu_k = \mu_k^0(T, P) + RT \ln a_k. \f]. - void getChemPotentials(doublereal* mu) const; - - /// Valid chemical potentials. Write into array \a mu the - /// chemical potentials of all species with thermo data valid - /// for the current temperature [J/kmol]. For other species, - /// set the chemical potential to the value \a not_mu. If \a - /// standard is set to true, then the values returned are - /// standard chemical potentials. - void getValidChemPotentials(doublereal not_mu, doublereal* mu, - bool standard = false) const; - - /// Temperature [K]. - doublereal temperature() const { return m_temp; } - - /// Set the mixture to a state of chemical equilibrium. - /// @param XY Integer flag specifying properties to hold fixed. - /// @param err Error tolerance for \f$\Delta \mu/RT \f$ for - /// all reactions. Also used as the relative error tolerance - /// for the outer loop. - /// @param maxsteps Maximum number of steps to take in solving - /// the fixed TP problem. - /// @param maxiter Maximum number of "outer" iterations for - /// problems holding fixed something other than (T,P). - /// @param loglevel Level of diagnostic output, written to a - /// file in HTML format. - doublereal equilibrate(int XY, doublereal err = 1.0e-9, - int maxsteps = 1000, int maxiter = 200, int loglevel = -99); - - - /// Set the temperature [K]. - void setTemperature(doublereal T) { - m_temp = T; - updatePhases(); - } - - /// Pressure [Pa]. - doublereal pressure() const { - return m_press; - } - - /// Volume [m^3]. - doublereal volume() const; - - /// Set the pressure [Pa]. - void setPressure(doublereal P) { - m_press = P; - updatePhases(); - } - - /// Enthalpy [J]. - doublereal enthalpy() const; - - /// Entropy [J/K]. - doublereal entropy() const; - - /// Gibbs function [J]. - doublereal gibbs() const; - - /// Heat capacity at constant pressure [J/K]. - doublereal cp() const; - - /// Number of phases. - index_t nPhases() const { - return m_np; - } - - /// Return true is species \a k is a species in a - /// multicomponent solution phase. - bool solutionSpecies(index_t k) const; - - index_t speciesPhaseIndex(index_t k) const{ - return m_spphase[k]; - } - - doublereal moleFraction(index_t k) const{ - return m_moleFractions[k]; - } - - void setPhaseMoleFractions(index_t n, doublereal* x); - - void setMolesByName(compositionMap& xMap); - - void setMolesByName(const std::string& x); - - void getMoles(doublereal * molNum) const; - void setMoles(doublereal* n); - - /// Return true if the phase \a p has valid thermo data for - /// the current temperature. - bool tempOK(index_t p) const { - return m_temp_OK[p]; - } - - protected: - - // These methods are meant for internal use. - - /// update the locally-stored composition to match the current - /// compositions of the phase objects. - void updateMoleFractions(); - - /// Set the states of the phase objects to the locally-stored - /// state. Note that if individual phases have T and P different - /// than that stored locally, the phase T and P will be modified. - void updatePhases() const; - - /** - * Vector of the number of moles in each phase. - * Length = m_np, number of phases. - */ - vector_fp m_moles; - - /** - * Vector of the ThermoPhase Pointers. - */ - std::vector m_phase; - array_t m_atoms; - /** - * Locally storred vector of mole fractions of all species - * comprising the MultiPhase object. - */ - vector_fp m_moleFractions; - vector_int m_spphase; - vector_int m_spstart; - std::vector m_enames; - vector_int m_atomicNumber; - std::vector m_snames; - mutable std::map m_enamemap; - /** - * Number of phases in the MultiPhase object - */ - index_t m_np; - doublereal m_temp; - doublereal m_press; - /** - * Number of distinct elements in all of the phases - */ - index_t m_nel; - /** - * Number of distinct species in all of the phases - */ - index_t m_nsp; - bool m_init; - int m_eloc; - mutable std::vector m_temp_OK; - doublereal m_Tmin, m_Tmax; - }; - - inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) { - size_t ip; - for (ip = 0; ip < x.nPhases(); ip++) { - if (x.phase(ip).name() != "") { - s << "*************** " << x.phase(ip).name() << " *****************" << std::endl; - } - else { - s << "*************** Phase " << ip << " *****************" << std::endl; - } - s << "Moles: " << x.phaseMoles(ip) << std::endl; - - s << report(x.phase(ip)) << std::endl; - } - return s; + /// Number of atoms of element \a m in species \a k. + doublereal nAtoms(int k, int m) { + if (!m_init) init(); + return m_atoms(m,k); } + /// Species mole fractions. Write the array of species mole + /// fractions into array \c x. The mole fractions are + /// normalized to sum to one in each phase. + void getMoleFractions(doublereal* x) const { + std::copy(m_moleFractions.begin(), m_moleFractions.end(), x); + } + + /// Process phases and build atomic composition array. After + /// init() has been called, no more phases may be added. + void init(); + + /// Moles of phase n. + doublereal phaseMoles(index_t n) const { + return m_moles[n]; + } + + /// Set the number of moles of phase with index n. + void setPhaseMoles(index_t n, doublereal moles) { + m_moles[n] = moles; + } + + /// Return a reference to phase n. The state of phase n is + /// also updated to match the state stored locally in the + /// mixture object. + phase_t& phase(index_t n); + + //! Moles of species \c k. + /*! + * Returns the moles of global species k. + * units = kmol + */ + doublereal speciesMoles(index_t k) const; + + /// Index of the species belonging to phase number \c p + /// with index \c k within the phase. + int speciesIndex(index_t k, index_t p) const { + return m_spstart[p] + k; + } + + /// Minimum temperature for which all solution phases have + /// valid thermo data. Stoichiometric phases are not + /// considered, since they may have thermo data only valid for + /// conditions for which they are stable. + doublereal minTemp() const { return m_Tmin; } + + /// Maximum temperature for which all solution phases have + /// valid thermo data. Stoichiometric phases are not + /// considered, since they may have thermo data only valid for + /// conditions for which they are stable. + doublereal maxTemp() const { return m_Tmax; } + + /// Total charge (Coulombs). + doublereal charge() const; + + /// Charge (Coulombs) of phase with index \a p. + doublereal phaseCharge(index_t p) const; + + /// Total moles of element \a m, summed over all phases. + doublereal elementMoles(index_t m) const; + + /// Chemical potentials. Write into array \a mu the chemical + /// potentials of all species [J/kmol]. The chemical + /// potentials are related to the activities by + /// \f[ \mu_k = \mu_k^0(T, P) + RT \ln a_k. \f]. + void getChemPotentials(doublereal* mu) const; + + /// Valid chemical potentials. Write into array \a mu the + /// chemical potentials of all species with thermo data valid + /// for the current temperature [J/kmol]. For other species, + /// set the chemical potential to the value \a not_mu. If \a + /// standard is set to true, then the values returned are + /// standard chemical potentials. + void getValidChemPotentials(doublereal not_mu, doublereal* mu, + bool standard = false) const; + + /// Temperature [K]. + doublereal temperature() const { return m_temp; } + + /// Set the mixture to a state of chemical equilibrium. + /// @param XY Integer flag specifying properties to hold fixed. + /// @param err Error tolerance for \f$\Delta \mu/RT \f$ for + /// all reactions. Also used as the relative error tolerance + /// for the outer loop. + /// @param maxsteps Maximum number of steps to take in solving + /// the fixed TP problem. + /// @param maxiter Maximum number of "outer" iterations for + /// problems holding fixed something other than (T,P). + /// @param loglevel Level of diagnostic output, written to a + /// file in HTML format. + doublereal equilibrate(int XY, doublereal err = 1.0e-9, + int maxsteps = 1000, int maxiter = 200, int loglevel = -99); + + + /// Set the temperature [K]. + void setTemperature(doublereal T) { + m_temp = T; + updatePhases(); + } + + /// Pressure [Pa]. + doublereal pressure() const { + return m_press; + } + + /// Volume [m^3]. + doublereal volume() const; + + /// Set the pressure [Pa]. + void setPressure(doublereal P) { + m_press = P; + updatePhases(); + } + + /// Enthalpy [J]. + doublereal enthalpy() const; + + /// Entropy [J/K]. + doublereal entropy() const; + + /// Gibbs function [J]. + doublereal gibbs() const; + + /// Heat capacity at constant pressure [J/K]. + doublereal cp() const; + + /// Number of phases. + index_t nPhases() const { + return m_np; + } + + /// Return true is species \a k is a species in a + /// multicomponent solution phase. + bool solutionSpecies(index_t k) const; + + //! Returns the phase index of the Kth "global" species + /*! + * @param k Global species index. + * + * @return + * Returns the index of the owning phase. + */ + index_t speciesPhaseIndex(index_t k) const { + return m_spphase[k]; + } + + //! Returns the mole fraction of global species k + doublereal moleFraction(index_t k) const{ + return m_moleFractions[k]; + } + + void setPhaseMoleFractions(index_t n, doublereal* x); + + void setMolesByName(compositionMap& xMap); + + void setMolesByName(const std::string& x); + + void getMoles(doublereal * molNum) const; + void setMoles(doublereal* n); + + //! Retrieves a vector of element abundances + /*! + * @param elemAbundances Vector of element abundances + * Length = number of elements in the MultiPhase object. + * Index is the global element index + * units is in kmol. + */ + void getElemAbundances(doublereal * elemAbundances) const; + + /// Return true if the phase \a p has valid thermo data for + /// the current temperature. + bool tempOK(index_t p) const { + return m_temp_OK[p]; + } + + protected: + + // These methods are meant for internal use. + + /// update the locally-stored composition to match the current + /// compositions of the phase objects. + void updateMoleFractions(); + + /// Set the states of the phase objects to the locally-stored + /// state. Note that if individual phases have T and P different + /// than that stored locally, the phase T and P will be modified. + void updatePhases() const; + + //! Calculate the element abundance vector + void calcElemAbundances() const; + /** + * Vector of the number of moles in each phase. + * Length = m_np, number of phases. + */ + vector_fp m_moles; + + /** + * Vector of the ThermoPhase Pointers. + */ + std::vector m_phase; + + //! Global Stoichiometric Coefficient array + /*! + * This is a two dimensional array m_atoms(m, k). The first + * index is the global element index. The second index, k, is the + * global species index. + * The value is the number of atoms of type m in species k. + */ + array_t m_atoms; + + /** + * Locally storred vector of mole fractions of all species + * comprising the MultiPhase object. + */ + vector_fp m_moleFractions; + vector_int m_spphase; + vector_int m_spstart; + std::vector m_enames; + vector_int m_atomicNumber; + + //! Vector of species names in the problem + /*! + * Vector is over all species defined in the object, + * the global species index. + */ + std::vector m_snames; + + mutable std::map m_enamemap; + /** + * Number of phases in the MultiPhase object + */ + index_t m_np; + doublereal m_temp; + doublereal m_press; + /** + * Number of distinct elements in all of the phases + */ + index_t m_nel; + /** + * Number of distinct species in all of the phases + */ + index_t m_nsp; + bool m_init; + int m_eloc; + mutable std::vector m_temp_OK; + doublereal m_Tmin, m_Tmax; + + mutable vector_fp m_elemAbundances; + }; + + inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) { + size_t ip; + for (ip = 0; ip < x.nPhases(); ip++) { + if (x.phase(ip).name() != "") { + s << "*************** " << x.phase(ip).name() << " *****************" << std::endl; + } + else { + s << "*************** Phase " << ip << " *****************" << std::endl; + } + s << "Moles: " << x.phaseMoles(ip) << std::endl; + + s << report(x.phase(ip)) << std::endl; + } + return s; + } + int BasisOptimize( int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, vector_int & orderVectorSpecies, diff --git a/tools/doc/Cantera.cfg.in b/tools/doc/Cantera.cfg.in index 95b3d318c..1bbf98571 100755 --- a/tools/doc/Cantera.cfg.in +++ b/tools/doc/Cantera.cfg.in @@ -122,7 +122,7 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \ WaterPropsIAPWS.h WaterPropsIAPWS.cpp \ WaterSSTP.h WaterSSTP.cpp \ PureFluidPhase.h PureFluidPhase.cpp \ - equil.h + equil.h MultiPhase.h MultiPhase.cpp RECURSIVE = NO EXCLUDE = CVS examples converters zeroD EXCLUDE_SYMLINKS = NO