diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h index 2bc73bad1..524e232de 100644 --- a/include/cantera/equil/MultiPhase.h +++ b/include/cantera/equil/MultiPhase.h @@ -5,7 +5,6 @@ */ // Copyright 2004 California Institute of Technology - #ifndef CT_MULTIPHASE_H #define CT_MULTIPHASE_H @@ -19,73 +18,55 @@ namespace Cantera //! A class for multiphase mixtures. The mixture can contain any //! number of phases of any type. /*! - * This object is the basic tool used by Cantera for use in - * Multiphase equilibrium calculations. + * This object is the basic tool used by Cantera for use in Multiphase + * equilibrium calculations. * - * It is a container for a set of phases. Each phase has a - * given number of kmoles. Therefore, MultiPhase may be considered - * an "extrinsic" thermodynamic object, in contrast to the ThermoPhase - * object, which is an "intrinsic" thermodynamic object. + * It is a container for a set of phases. Each phase has a given number of + * kmoles. Therefore, MultiPhase may be considered an "extrinsic" + * thermodynamic object, in contrast to the ThermoPhase object, which is an + * "intrinsic" thermodynamic object. * - * MultiPhase may be considered to be "upstream" of the ThermoPhase - * objects in the sense that setting a property within MultiPhase, - * such as temperature, pressure, or species mole number, - * affects the underlying ThermoPhase object, but not the - * other way around. + * MultiPhase may be considered to be "upstream" of the ThermoPhase objects in + * the sense that setting a property within MultiPhase, such as temperature, + * pressure, or species mole number, affects the underlying ThermoPhase + * object, but not the other way around. * - * All phases have the same - * temperature and pressure, and a specified number of moles for - * each phase. - * 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 intersection of the elements of each phase. + * All phases have the same temperature and pressure, and a specified number + * of moles for each phase. 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 intersection of + * the elements of each phase. * - * Below, reference is made to global species and global elements. - * These refer to the collective species and elements encompassing - * all of the phases tracked by the object. + * Below, reference is made to global species and global elements. These refer + * to the collective species and elements encompassing all of the phases + * tracked by the object. * - * The global element list kept by this object is an - * intersection of the element lists of all the phases that - * comprise the MultiPhase. + * The global element list kept by this object is an intersection of the + * element lists of all the phases that comprise the MultiPhase. * - * The global species list kept by this object is a - * concatenated list of all of the species in all the phases that - * comprise the MultiPhase. The ordering of species is contiguous - * with respect to the phase id. + * The global species list kept by this object is a concatenated list of all + * of the species in all the phases that comprise the MultiPhase. The ordering + * of species is contiguous with respect to the phase id. * - * @ingroup equilfunctions + * @ingroup equilfunctions */ class MultiPhase { - public: //! Constructor. /*! - * The constructor takes no arguments, since - * phases are added using method addPhase(). + * The constructor takes no arguments, since phases are added using + * method addPhase(). */ MultiPhase(); - //! Copy Constructor - /*! - * @param right Object to be copied - */ MultiPhase(const MultiPhase& right); - //! Destructor. - /*! - * Does nothing. Class MultiPhase does not take - * "ownership" (i.e. responsibility for destroying) the - * phase objects. - */ + //! Destructor. Does nothing. Class MultiPhase does not take "ownership" + //! (i.e. responsibility for destroying) the phase objects. virtual ~MultiPhase(); - //! Assignment operator - /*! - * @param right Object to be copied - */ MultiPhase& operator=(const MultiPhase& right); //! Add a vector of phases to the mixture @@ -114,21 +95,21 @@ public: */ void addPhase(ThermoPhase* p, doublereal moles); - /// Number of elements. + //! Number of elements. size_t nElements() const { return m_nel; } - //! Check that the specified element index is in range + //! Check that the specified element index is in range. //! Throws an exception if m is greater than nElements()-1 void checkElementIndex(size_t m) const; - //! Check that an array size is at least nElements() + //! Check that an array size is at least nElements(). //! Throws an exception if mm is less than nElements(). Used before calls //! which take an array pointer. void checkElementArraySize(size_t mm) const; - //! Returns the string name of the global element \a m. + //! Returns the name of the global element *m*. /*! * @param m index of the global element */ @@ -145,11 +126,11 @@ public: return m_nsp; } - //! Check that the specified species index is in range + //! Check that the specified species index is in range. //! Throws an exception if k is greater than nSpecies()-1 void checkSpeciesIndex(size_t k) const; - //! Check that an array size is at least nSpecies() + //! Check that an array size is at least nSpecies(). //! Throws an exception if kk is less than nSpecies(). Used before calls //! which take an array pointer. void checkSpeciesArraySize(size_t kk) const; @@ -181,10 +162,10 @@ public: void getMoleFractions(doublereal* const x) const; //! Process phases and build atomic composition array. - /*!This method - * must be called after all phases are added, before doing - * anything else with the mixture. After init() has been called, - * no more phases may be added. + /*! + * This method must be called after all phases are added, before doing + * anything else with the mixture. After init() has been called, no more + * phases may be added. */ void init(); @@ -216,14 +197,13 @@ public: */ void setPhaseMoles(const size_t n, const doublereal moles); - /// Return a %ThermoPhase reference to phase n. - /*! The state of phase n is - * also updated to match the state stored locally in the - * mixture object. + /// Return a reference to phase n. + /*! + * The state of phase n is also updated to match the state stored locally + * in the mixture object. * * @param n Phase Index - * - * @return Reference to the %ThermoPhase object for the phase + * @return Reference to the ThermoPhase object for the phase */ thermo_t& phase(size_t n); @@ -236,11 +216,8 @@ public: //! which take an array pointer. void checkPhaseArraySize(size_t mm) const; - //! Returns the moles of global species \c k. + //! Returns the moles of global species \c k. units = kmol /*! - * Returns the moles of global species k. - * units = kmol - * * @param kGlob Global species index k */ doublereal speciesMoles(size_t kGlob) const; @@ -248,8 +225,6 @@ public: //! Return the global index of the species belonging to phase number \c p //! with local index \c k within the phase. /*! - * Returns the index of the global species - * * @param k local index of the species within the phase * @param p index of the phase */ @@ -260,8 +235,6 @@ public: //! Return the global index of the species belonging to phase name \c phaseName //! with species name \c speciesName /*! - * Returns the index of the global species - * * @param speciesName Species Name * @param phaseName Phase Name * @@ -272,18 +245,16 @@ public: */ size_t speciesIndex(const std::string& speciesName, const std::string& phaseName); - /// 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. + /// 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. + /// 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; } @@ -293,7 +264,9 @@ public: /// Charge (Coulombs) of phase with index \a p. /*! - * @param p Phase Index + * The net charge is computed as \f[ Q_p = N_p \sum_k F z_k X_k \f] + * where the sum runs only over species in phase \a p. + * @param p index of the phase for which the charge is desired. */ doublereal phaseCharge(size_t p) const; @@ -305,61 +278,55 @@ public: //! Returns a vector of Chemical potentials. /*! - * Write into array \a mu the chemical - * potentials of all species [J/kmol]. The chemical - * potentials are related to the activities by + * 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$. * - * @param mu Chemical potential vector. - * Length = num global species. - * Units = J/kmol. + * @param mu Chemical potential vector. Length = num global species. Units + * = J/kmol. */ void getChemPotentials(doublereal* mu) const; /// Returns a vector of 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. + * 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. * - * This method is designed for use in computing chemical - * equilibrium by Gibbs minimization. For solution phases (more - * than one species), this does the same thing as - * getChemPotentials. But for stoichiometric phases, this writes - * into array \a mu the user-specified value \a not_mu instead of - * the chemical potential if the temperature is outside the range - * for which the thermo data for the one species in the phase are - * valid. The need for this arises since many condensed phases - * have thermo data fit only for the temperature range for which - * they are stable. For example, in the NASA database, the fits - * for H2O(s) are only done up to 0 C, the fits for H2O(L) are - * only done from 0 C to 100 C, etc. Using the polynomial fits outside - * the range for which the fits were done can result in spurious - * chemical potentials, and can lead to condensed phases - * appearing when in fact they should be absent. + * This method is designed for use in computing chemical equilibrium by + * Gibbs minimization. For solution phases (more than one species), this + * does the same thing as getChemPotentials. But for stoichiometric + * phases, this writes into array \a mu the user-specified value \a not_mu + * instead of the chemical potential if the temperature is outside the + * range for which the thermo data for the one species in the phase are + * valid. The need for this arises since many condensed phases have thermo + * data fit only for the temperature range for which they are stable. For + * example, in the NASA database, the fits for H2O(s) are only done up to + * 0 C, the fits for H2O(L) are only done from 0 C to 100 C, etc. Using + * the polynomial fits outside the range for which the fits were done can + * result in spurious chemical potentials, and can lead to condensed + * phases appearing when in fact they should be absent. * - * By setting \a not_mu to a large positive value, it is possible - * to force routines which seek to minimize the Gibbs free energy - * of the mixture to zero out any phases outside the temperature - * range for which their thermo data are valid. + * By setting \a not_mu to a large positive value, it is possible to force + * routines which seek to minimize the Gibbs free energy of the mixture to + * zero out any phases outside the temperature range for which their + * thermo data are valid. * - * @param not_mu Value of the chemical potential to set - * species in phases, for which the thermo data - * is not valid + * @param not_mu Value of the chemical potential to set species in phases, + * for which the thermo data is not valid * - * @param mu Vector of chemical potentials - * length = Global species, units = J kmol-1 + * @param mu Vector of chemical potentials. length = Global species, + * units = J kmol-1 * - * @param standard If this method is called with \a standard set to true, then - * the composition-independent standard chemical potentials are - * returned instead of the composition-dependent chemical - * potentials. + * @param standard If this method is called with \a standard set to true, + * then the composition-independent standard chemical + * potentials are returned instead of the composition- + * dependent chemical potentials. */ void getValidChemPotentials(doublereal not_mu, doublereal* mu, bool standard = false) const; @@ -371,21 +338,20 @@ public: //! 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. + * @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]. /*! * @param T value of the temperature (Kelvin) @@ -396,7 +362,6 @@ public: /*! * @param T Temperature of the system (kelvin) * @param Pres pressure of the system (pascal) - * (kmol) */ void setState_TP(const doublereal T, const doublereal Pres); @@ -414,38 +379,39 @@ public: return m_press; } - /// Volume [m^3]. + /// The total mixture volume [m^3]. /*! * Returns the cumulative sum of the volumes of all the - * phases in the %MultiPhase. + * phases in the mixture. */ doublereal volume() const; //! Set the pressure [Pa]. /*! - * @param P Set the pressure in the %MultiPhase object (Pa) + * @param P Set the pressure in the MultiPhase object (Pa) */ void setPressure(doublereal P) { m_press = P; updatePhases(); } - /// Enthalpy [J]. + //! The enthalpy of the mixture [J]. doublereal enthalpy() const; - /// Enthalpy [J]. + //! The internal energy of the mixture [J]. doublereal IntEnergy() const; - /// Entropy [J/K]. + //! The entropy of the mixture [J/K]. doublereal entropy() const; - /// Gibbs function [J]. + //! The Gibbs function of the mixture [J]. doublereal gibbs() const; - /// Heat capacity at constant pressure [J/K]. + //! Heat capacity at constant pressure [J/K]. Note that this does not + //! account for changes in composition of the mixture with temperature. doublereal cp() const; - /// Number of phases. + //! Number of phases. size_t nPhases() const { return m_np; } @@ -461,8 +427,7 @@ public: /*! * @param kGlob Global species index. * - * @return - * Returns the index of the owning phase. + * @return Returns the index of the owning phase. */ size_t speciesPhaseIndex(const size_t kGlob) const; @@ -474,58 +439,49 @@ public: //! Set the Mole fractions of the nth phase /*! - * This function sets the mole fractions of the - * nth phase. Note, the mole number of the phase - * stays constant + * This function sets the mole fractions of the nth phase. Note, the mole + * number of the phase stays constant * - * @param n ID of the phase + * @param n index of the phase * @param x Vector of input mole fractions. */ void setPhaseMoleFractions(const size_t n, const doublereal* const x); - //! Set the number numbers of species in the MultiPhase + //! Set the number of moles of species in the mixture /*! - * @param xMap CompositionMap of the species with - * nonzero mole numbers - * units = kmol. + * @param xMap CompositionMap of the species with nonzero mole numbers. + * Mole numbers that are less than or equal to zero will be + * set to zero. units = kmol. */ void setMolesByName(compositionMap& xMap); - //! Set the Moles via a string containing their names. + //! Set the moles via a string containing their names. /*! - * The string x is in the form of a composition map - * Species which are not listed by name in the composition - * map are set to zero. + * The string x is in the form of a composition map. Species which are not + * listed are set to zero. * * @param x string x in the form of a composition map * where values are the moles of the species. */ void setMolesByName(const std::string& x); - - //! Return a vector of global species mole numbers + //! Get the mole numbers of all species in the multiphase object /*! - * Returns a vector of the number of moles of each species - * in the multiphase object. - * - * @param molNum Vector of doubles of length nSpecies - * containing the global mole numbers - * (kmol). + * @param[out] molNum Vector of doubles of length nSpecies containing the + * global mole numbers (kmol). */ void getMoles(doublereal* molNum) const; //! Sets all of the global species mole numbers /*! - * Sets the number of moles of each species - * in the multiphase object. + * The state of each phase object is also updated to have the specified + * composition and the mixture temperature and pressure. * - * @param n Vector of doubles of length nSpecies - * containing the global mole numbers - * (kmol). + * @param n Vector of doubles of length nSpecies containing the global + * mole numbers (kmol). */ void setMoles(const doublereal* n); - //! Adds moles of a certain species to the mixture /*! * @param indexS Index of the species in the MultiPhase object @@ -537,8 +493,7 @@ public: /*! * @param elemAbundances Vector of element abundances * Length = number of elements in the MultiPhase object. - * Index is the global element index - * units is in kmol. + * Index is the global element index. Units is in kmol. */ void getElemAbundances(doublereal* elemAbundances) const; @@ -554,30 +509,28 @@ public: //! Update the locally-stored composition within this object //! to match the current compositions of the phase objects. /*! - * Query the underlying ThermoPhase objects for their mole - * fractions and fill in the mole fraction vector of this - * current object. Adjust element compositions within this - * object to match. + * Query the underlying ThermoPhase objects for their mole fractions and + * fill in the mole fraction vector of this current object. Adjust + * element compositions within this object to match. * - * This is an upload operation in the sense that we are taking - * downstream information (ThermoPhase object info) and - * applying it to an upstream object (MultiPhase object). + * This is an upload operation in the sense that we are taking downstream + * information (ThermoPhase object info) and applying it to an upstream + * object (MultiPhase object). */ void uploadMoleFractionsFromPhases(); //! Set the states of the phase objects to the locally-stored //! state within this MultiPhase object. /*! + * This method sets each phase to the mixture temperature and pressure, + * and sets the phase mole fractions based on the mixture mole numbers. * - * Note that if individual phases have T and P different - * than that stored locally, the phase T and P will be modified. + * This is an download operation in the sense that we are taking + * upstream object information (MultiPhase object) and + * applying it to downstream objects (ThermoPhase object information) * - * This is an download operation in the sense that we are taking - * upstream object information (MultiPhase object) and - * applying it to downstream objects (ThermoPhase object information) - * - * Therefore, the term, "update", is appropriate for a downstream - * operation. + * Therefore, the term, "update", is appropriate for a downstream + * operation. */ void updatePhases() const; @@ -585,24 +538,20 @@ private: //! 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. - */ + //! 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. + * 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. */ DenseMatrix m_atoms; @@ -619,11 +568,11 @@ private: */ std::vector m_spphase; - //! Vector of ints containing of first species index in the global list of species - //! for each phase + //! Vector of ints containing of first species index in the global list of + //! species for each phase /*! - * kfirst = m_spstart[ip], kfirst is the index of the first species in the ip'th - * phase. + * kfirst = m_spstart[ip], kfirst is the index of the first species in + * the ip'th phase. */ std::vector m_spstart; @@ -652,9 +601,7 @@ private: */ std::map m_enamemap; - /** - * Number of phases in the MultiPhase object - */ + //! Number of phases in the MultiPhase object size_t m_np; //! Current value of the temperature (kelvin) @@ -663,13 +610,10 @@ private: //! Current value of the pressure (Pa) doublereal m_press; - /** - * Number of distinct elements in all of the phases - */ + //! Number of distinct elements in all of the phases size_t m_nel; - /** - * Number of distinct species in all of the phases - */ + + //! Number of distinct species in all of the phases size_t m_nsp; //! True if the init() routine has been called, and the MultiPhase frozen @@ -697,8 +641,7 @@ private: //! Minimum temperature for which thermo parameterizations are valid /*! - * Stoichiometric phases are ignored in this determination. - * units Kelvin + * Stoichiometric phases are ignored in this determination. units Kelvin */ doublereal m_Tmax; @@ -735,13 +678,11 @@ inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) return s; } -//! Choose the optimum basis of species for the equilibrium calculations. +//! Choose the optimum basis of species for the equilibrium calculations. /*! - * This is done by - * choosing the species with the largest mole fraction - * not currently a linear combination of the previous components. - * Then, calculate the stoichiometric coefficient matrix for that - * basis. + * This is done by choosing the species with the largest mole fraction not + * currently a linear combination of the previous components. Then, calculate + * the stoichiometric coefficient matrix for that basis. * * Calculates the identity of the component species in the mechanism. * Rearranges the solution data to put the component data at the @@ -750,87 +691,65 @@ inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) * Then, calculates SC(J,I) the formation reactions for all noncomponent * species in the mechanism. * - * Input - * --------- - * @param mphase Pointer to the multiphase object. Contains the - * species mole fractions, which are used to pick the - * current optimal species component basis. - * @param orderVectorElements - * Order vector for the elements. The element rows - * in the formula matrix are - * rearranged according to this vector. - * @param orderVectorSpecies - * Order vector for the species. The species are - * rearranged according to this formula. The first - * nCompoments of this vector contain the calculated - * species components on exit. - * @param doFormRxn If true, the routine calculates the formation - * reaction matrix based on the calculated - * component species. If false, this step is skipped. + * @param[in] mphase Pointer to the multiphase object. Contains the species + * mole fractions, which are used to pick the current optimal species + * component basis. + * @param[in] orderVectorElements Order vector for the elements. The element + * rows in the formula matrix are rearranged according to this vector. + * @param[in] orderVectorSpecies Order vector for the species. The species are + * rearranged according to this formula. The first nCompoments of this + * vector contain the calculated species components on exit. + * @param[in] doFormRxn If true, the routine calculates the formation + * reaction matrix based on the calculated component species. If + * false, this step is skipped. + * @param[out] usedZeroedSpecies = If true, then a species with a zero + * concentration was used as a component. The problem may be + * converged. + * @param[out] formRxnMatrix + * @return The number of components. * - * Output - * --------- - * @param usedZeroedSpecies = If true, then a species with a zero concentration - * was used as a component. The problem may be - * converged. - * @param formRxnMatrix - * - * @return Returns the number of components. - * - * @ingroup equilfunctions + * @ingroup equilfunctions */ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, std::vector& orderVectorSpecies, std::vector& orderVectorElements, vector_fp& formRxnMatrix); -//! This subroutine handles the potential rearrangement of the constraint -//! equations represented by the Formula Matrix. +//! Handles the potential rearrangement of the constraint equations +//! represented by the Formula Matrix. /*! - * Rearrangement is only - * necessary when the number of components is less than the number of - * elements. For this case, some constraints can never be satisfied - * exactly, because the range space represented by the Formula - * Matrix of the components can't span the extra space. These - * constraints, which are out of the range space of the component - * Formula matrix entries, are migrated to the back of the Formula - * matrix. + * Rearrangement is only necessary when the number of components is less + * than the number of elements. For this case, some constraints can never + * be satisfied exactly, because the range space represented by the Formula + * Matrix of the components can't span the extra space. These constraints, + * which are out of the range space of the component Formula matrix + * entries, are migrated to the back of the Formula matrix. * - * A prototypical example is an extra element column in - * FormulaMatrix[], - * which is identically zero. For example, let's say that argon is - * has an element column in FormulaMatrix[], but no species in the - * mechanism - * actually contains argon. Then, nc < ne. Unless the entry for - * desired element abundance vector for Ar is zero, then this - * element abundance constraint can never be satisfied. The - * constraint vector is not in the range space of the formula - * matrix. - * Also, without perturbation - * of FormulaMatrix[], BasisOptimize[] would produce a zero pivot - * because the matrix - * would be singular (unless the argon element column was already the - * last column of FormulaMatrix[]. - * This routine borrows heavily from BasisOptimize algorithm. It - * finds nc constraints which span the range space of the Component - * Formula matrix, and assigns them as the first nc components in the - * formula matrix. This guarantees that BasisOptimize has a - * nonsingular matrix to invert. - * input - * @param nComponents Number of components calculated previously. + * A prototypical example is an extra element column in FormulaMatrix[], which + * is identically zero. For example, let's say that argon is has an element + * column in FormulaMatrix[], but no species in the mechanism actually + * contains argon. Then, nc < ne. Unless the entry for desired element + * abundance vector for Ar is zero, then this element abundance constraint can + * never be satisfied. The constraint vector is not in the range space of the + * formula matrix. * - * @param elementAbundances Current value of the element abundances + * Also, without perturbation of FormulaMatrix[], BasisOptimize[] would + * produce a zero pivot because the matrix would be singular (unless the argon + * element column was already the last column of FormulaMatrix[]. * - * @param mphase Input pointer to a MultiPhase object + * This routine borrows heavily from BasisOptimize algorithm. It finds nc + * constraints which span the range space of the Component Formula matrix, and + * assigns them as the first nc components in the formula matrix. This + * guarantees that BasisOptimize has a nonsingular matrix to invert. * - * @param orderVectorSpecies input vector containing the ordering - * of the global species in mphase. This is used - * to extract the component basis of the mphase object. - * - * output - * @param orderVectorElements Output vector containing the order - * of the elements that is necessary for - * calculation of the formula matrix. + * @param[in] nComponents Number of components calculated previously. + * @param[in] elementAbundances Current value of the element abundances + * @param[in] mphase Input pointer to a MultiPhase object + * @param[in] orderVectorSpecies input vector containing the ordering of the + * global species in mphase. This is used to extract the component + * basis of the mphase object. + * @param[out] orderVectorElements Output vector containing the order of the + * elements that is necessary for calculation of the formula matrix. * * @ingroup equilfunctions */ diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index a3ba6fc5d..c8ecbf7f1 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -1,6 +1,5 @@ /** - * @file BasisOptimize.cpp - * Functions which calculation optimized basis of the + * @file BasisOptimize.cpp Functions which calculation optimized basis of the * stoichiometric coefficient matrix (see /ref equil functions) */ #include "cantera/base/ct_defs.h" @@ -19,10 +18,10 @@ namespace Cantera int BasisOptimize_print_lvl = 0; } -//! Print a string within a given space limit. This routine limits the amount of the string that will be printed to a -//! maximum of "space" characters. +//! Print a string within a given space limit. /*! - * + * This routine limits the amount of the string that will be printed to a + * maximum of "space" characters. * @param str String -> must be null terminated. * @param space space limit for the printing. * @param alignment 0 centered @@ -32,98 +31,46 @@ int BasisOptimize_print_lvl = 0; static void print_stringTrunc(const char* str, int space, int alignment); #endif - -//! Finds the location of the maximum component in a double vector INPUT +//! Finds the location of the maximum component in a vector *x* /*! - * @param x Vector to search - * @param j j <= i < n : i is the range of indices to search in X(*) - * @param n Length of the vector + * @param x Vector to search + * @param j j <= i < n : i is the range of indices to search in *x* + * @param n Length of the vector * - * @return index of the greatest value on X(*) searched + * @return index of the greatest value on *x* searched */ static size_t amax(double* x, size_t j, size_t n); -//! Invert an nxn matrix and solve m rhs's +//! Invert an nxn matrix and solve m rhs's /*! + * Solve C X + B = 0 * - * Solve C X + B = 0; + * This routine uses Gauss elimination and is optimized for the solution of + * lots of rhs's. A crude form of row pivoting is used here. * - * This routine uses Gauss elimination and is optimized for the solution - * of lots of rhs's. - * A crude form of row pivoting is used here. + * @param c C is the matrix to be inverted + * @param idem first dimension in the calling routine. + * idem >= n must be true + * @param n number of rows and columns in the matrix + * @param b rhs of the matrix problem + * @param m number of rhs to be solved for * - * @param c C is the matrix to be inverted - * @param idem first dimension in the calling routine - * idem >= n must be true - * @param n number of rows and columns in the matrix - * @param b rhs of the matrix problem - * @param m number of rhs to be solved for + * - c[i+j*idem] = c_i_j = Matrix to be inverted + * - b[i+j*idem] = b_i_j = vectors of rhs's. Each column is a new rhs. * - * c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number - * j = column number - * b[i+j*idem] = b_i_j = vectors of rhs's: i = row number - * j = column number - * (each column is a new rhs) + * Where j = column number and i = row number. * - * @return Retuns the value - * 1 : Matrix is singular - * 0 : solution is OK + * @return Retuns 1 if the matrix is singular, or 0 if the solution is OK * - * The solution is returned in the matrix b. + * The solution is returned in the matrix b. */ static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m); -/* - * Choose the optimum basis for the calculations. This is done by - * choosing the species with the largest mole fraction - * not currently a linear combination of the previous components. - * Then, calculate the stoichiometric coefficient matrix for that - * basis. - * - * Calculates the identity of the component species in the mechanism. - * Rearranges the solution data to put the component data at the - * front of the species list. - * - * Then, calculates SC(J,I) the formation reactions for all noncomponent - * species in the mechanism. - * - * Input - * --------- - * mphase Pointer to the multiphase object. Contains the - * species mole fractions, which are used to pick the - * current optimal species component basis. - * orderVectorElement - * Order vector for the elements. The element rows - * in the formula matrix are - * rearranged according to this vector. - * orderVectorSpecies - * Order vector for the species. The species are - * rearranged according to this formula. The first - * nCompoments of this vector contain the calculated - * species components on exit. - * doFormRxn If true, the routine calculates the formation - * reaction matrix based on the calculated - * component species. If false, this step is skipped. - * - * Output - * --------- - * usedZeroedSpecies = If true, then a species with a zero concentration - * was used as a component. The problem may be - * converged. - * formRxnMatrix - * - * Return - * -------------- - * returns the number of components. - * - * - */ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, std::vector& orderVectorSpecies, std::vector& orderVectorElements, vector_fp& formRxnMatrix) { - size_t j, jj, k=0, kk, l, i, jl, ml; bool lindep; std::string ename; @@ -547,34 +494,6 @@ static size_t amax(double* x, size_t j, size_t n) return largest; } -/* - * vcs_mlequ: - * - * Invert an nxn matrix and solve m rhs's - * - * Solve C X + B = 0; - * - * This routine uses Gauss elimination and is optimized for the solution - * of lots of rhs's. - * A crude form of row pivoting is used here. - * - * - * c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number - * j = column number - * b[i+j*idem] = b_i_j = vectors of rhs's: i = row number - * j = column number - * (each column is a new rhs) - * n = number of rows and columns in the matrix - * m = number of rhs to be solved for - * idem = first dimension in the calling routine - * idem >= n must be true - * - * Return Value - * 1 : Matrix is singular - * 0 : solution is OK - * - * The solution is returned in the matrix b. - */ static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m) { size_t i, j, k, l; @@ -637,50 +556,13 @@ static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m) } } return 0; -} /* mlequ() *************************************************************/ +} - -/* - * - * ElemRearrange: - * - * This subroutine handles the rearrangement of the constraint - * equations represented by the Formula Matrix. Rearrangement is only - * necessary when the number of components is less than the number of - * elements. For this case, some constraints can never be satisfied - * exactly, because the range space represented by the Formula - * Matrix of the components can't span the extra space. These - * constraints, which are out of the range space of the component - * Formula matrix entries, are migrated to the back of the Formula - * matrix. - * - * A prototypical example is an extra element column in - * FormulaMatrix[], - * which is identically zero. For example, let's say that argon is - * has an element column in FormulaMatrix[], but no species in the - * mechanism - * actually contains argon. Then, nc < ne. Unless the entry for - * desired element abundance vector for Ar is zero, then this - * element abundance constraint can never be satisfied. The - * constraint vector is not in the range space of the formula - * matrix. - * Also, without perturbation - * of FormulaMatrix[], BasisOptimize[] would produce a zero pivot - * because the matrix - * would be singular (unless the argon element column was already the - * last column of FormulaMatrix[]. - * This routine borrows heavily from BasisOptimize algorithm. It - * finds nc constraints which span the range space of the Component - * Formula matrix, and assigns them as the first nc components in the - * formula matrix. This guarantees that BasisOptimize has a - * nonsingular matrix to invert. - */ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, MultiPhase* mphase, std::vector& orderVectorSpecies, std::vector& orderVectorElements) { - size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0; bool lindep = false; @@ -892,4 +774,4 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda */ } while (jr < (nComponents-1)); return nComponents; -} /* vcs_elem_rearrange() ****************************************************/ +} diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index bb318fa20..7013eeccd 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -16,8 +16,6 @@ using namespace std; namespace Cantera { -//==================================================================================================================== -// Constructor. MultiPhase::MultiPhase() : m_np(0), m_temp(298.15), @@ -30,11 +28,7 @@ MultiPhase::MultiPhase() : m_Tmax(100000.0) { } -//==================================================================================================================== -// Copy Constructor -/* - * @param right Object to be copied - */ + MultiPhase::MultiPhase(const MultiPhase& right) : m_np(0), m_temp(298.15), @@ -48,21 +42,11 @@ MultiPhase::MultiPhase(const MultiPhase& right) : { operator=(right); } -//==================================================================================================================== -// Destructor. -/* - * Does nothing. Class MultiPhase does not take - * "ownership" (i.e. responsibility for destroying) the - * phase objects. - */ + MultiPhase::~MultiPhase() { } -//==================================================================================================================== -// Assignment operator -/* - * @param right Object to be copied - */ + MultiPhase& MultiPhase::operator=(const MultiPhase& right) { if (&right != this) { @@ -89,7 +73,7 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right) } return *this; } -//==================================================================================================================== + void MultiPhase:: addPhases(MultiPhase& mix) { @@ -98,7 +82,7 @@ addPhases(MultiPhase& mix) addPhase(mix.m_phase[n], mix.m_moles[n]); } } -//==================================================================================================================== + void MultiPhase:: addPhases(std::vector& phases, const vector_fp& phaseMoles) { @@ -109,7 +93,7 @@ addPhases(std::vector& phases, const vector_fp& phaseMoles) } init(); } -//==================================================================================================================== + void MultiPhase:: addPhase(ThermoPhase* p, doublereal moles) { @@ -185,11 +169,7 @@ addPhase(ThermoPhase* p, doublereal moles) } } } -//==================================================================================================================== -// Process phases and build atomic composition array. This method -// must be called after all phases are added, before doing -// anything else with the mixture. After init() has been called, -// no more phases may be added. + void MultiPhase::init() { if (m_init) { @@ -253,10 +233,6 @@ void MultiPhase::init() updatePhases(); } -//==================================================================================================================== -// Return a reference to phase n. The state of phase n is -// also updated to match the state stored locally in the -// mixture object. ThermoPhase& MultiPhase::phase(size_t n) { if (!m_init) { @@ -282,18 +258,12 @@ void MultiPhase::checkPhaseArraySize(size_t mm) const } } -//==================================================================================================================== -/// Moles of species \c k. doublereal MultiPhase::speciesMoles(size_t k) const { size_t ip = m_spphase[k]; return m_moles[ip]*m_moleFractions[k]; } -//==================================================================================================================== -// Total moles of global element \a m, summed over all phases. -/* - * @param m Index of the global element - */ + doublereal MultiPhase::elementMoles(size_t m) const { doublereal sum = 0.0, phasesum; @@ -309,8 +279,7 @@ doublereal MultiPhase::elementMoles(size_t m) const } return sum; } -//==================================================================================================================== -// Total charge, summed over all phases + doublereal MultiPhase::charge() const { doublereal sum = 0.0; @@ -320,7 +289,7 @@ doublereal MultiPhase::charge() const } return sum; } -//==================================================================================================================== + size_t MultiPhase::speciesIndex(const std::string& speciesName, const std::string& phaseName) { if (!m_init) { @@ -336,11 +305,7 @@ size_t MultiPhase::speciesIndex(const std::string& speciesName, const std::strin } return m_spstart[p] + k; } -//==================================================================================================================== -/// Net charge of one phase (Coulombs). The net charge is computed as -/// \f[ Q_p = N_p \sum_k F z_k X_k \f] -/// where the sum runs only over species in phase \a p. -/// @param p index of the phase for which the charge is desired. + doublereal MultiPhase::phaseCharge(size_t p) const { doublereal phasesum = 0.0; @@ -351,9 +316,7 @@ doublereal MultiPhase::phaseCharge(size_t p) const } return Faraday*phasesum*m_moles[p]; } -//==================================================================================================================== -/// Get the chemical potentials of all species in all phases. void MultiPhase::getChemPotentials(doublereal* mu) const { size_t i, loc = 0; @@ -363,34 +326,7 @@ void MultiPhase::getChemPotentials(doublereal* mu) const loc += m_phase[i]->nSpecies(); } } -//==================================================================================================================== -// Get chemical potentials of species with valid thermo -// data. This method is designed for use in computing chemical -// equilibrium by Gibbs minimization. For solution phases (more -// than one species), this does the same thing as -// getChemPotentials. But for stoichiometric phases, this writes -// into array \a mu the user-specified value \a not_mu instead of -// the chemical potential if the temperature is outside the range -// for which the thermo data for the one species in the phase are -// valid. The need for this arises since many condensed phases -// have thermo data fit only for the temperature range for which -// they are stable. For example, in the NASA database, the fits -// for H2O(s) are only done up to 0 C, the fits for H2O(L) are -// only done from 0 C to 100 C, etc. Using the polynomial fits outside -// the range for which the fits were done can result in spurious -// chemical potentials, and can lead to condensed phases -// appearing when in fact they should be absent. -// -// By setting \a not_mu to a large positive value, it is possible -// to force routines which seek to minimize the Gibbs free energy -// of the mixture to zero out any phases outside the temperature -// range for which their thermo data are valid. -// -// If this method is called with \a standard set to true, then -// the composition-independent standard chemical potentials are -// returned instead of the composition-dependent chemical -// potentials. -// + void MultiPhase::getValidChemPotentials(doublereal not_mu, doublereal* mu, bool standard) const { @@ -411,8 +347,7 @@ void MultiPhase::getValidChemPotentials(doublereal not_mu, loc += m_phase[i]->nSpecies(); } } -//==================================================================================================================== -/// True if species \a k belongs to a solution phase. + bool MultiPhase::solutionSpecies(size_t k) const { if (m_phase[m_spphase[k]]->nSpecies() > 1) { @@ -421,8 +356,7 @@ bool MultiPhase::solutionSpecies(size_t k) const return false; } } -//==================================================================================================================== -/// The Gibbs free energy of the mixture (J). + doublereal MultiPhase::gibbs() const { size_t i; @@ -435,8 +369,7 @@ doublereal MultiPhase::gibbs() const } return sum; } -//==================================================================================================================== -/// The enthalpy of the mixture (J). + doublereal MultiPhase::enthalpy() const { size_t i; @@ -449,8 +382,7 @@ doublereal MultiPhase::enthalpy() const } return sum; } -//==================================================================================================================== -/// The internal energy of the mixture (J). + doublereal MultiPhase::IntEnergy() const { size_t i; @@ -463,8 +395,7 @@ doublereal MultiPhase::IntEnergy() const } return sum; } -//==================================================================================================================== -/// The entropy of the mixture (J/K). + doublereal MultiPhase::entropy() const { size_t i; @@ -477,10 +408,7 @@ doublereal MultiPhase::entropy() const } return sum; } -//==================================================================================================================== -/// The specific heat at constant pressure and composition (J/K). -/// Note that this does not account for changes in composition of -/// the mixture with temperature. + doublereal MultiPhase::cp() const { size_t i; @@ -494,10 +422,6 @@ doublereal MultiPhase::cp() const return sum; } -//==================================================================================================================== - -/// Set the mole fractions of phase \a n to the values in -/// array \a x. void MultiPhase::setPhaseMoleFractions(const size_t n, const doublereal* const x) { if (!m_init) { @@ -510,10 +434,7 @@ void MultiPhase::setPhaseMoleFractions(const size_t n, const doublereal* const x m_moleFractions[istart+k] = x[k]; } } -//==================================================================================================================== -// Set the species moles using a map. The map \a xMap maps -// species name strings to mole numbers. Mole numbers that are -// less than or equal to zero will be set to zero. + void MultiPhase::setMolesByName(compositionMap& xMap) { size_t kk = nSpecies(); @@ -527,18 +448,14 @@ void MultiPhase::setMolesByName(compositionMap& xMap) } setMoles(DATA_PTR(moles)); } -//==================================================================================================================== -// Set the species moles using a string. Unspecified species are -// set to zero. + void MultiPhase::setMolesByName(const std::string& x) { // build the composition map from the string, and then set the moles. compositionMap xx = parseCompString(x, m_snames); setMolesByName(xx); } -//==================================================================================================================== -// Get the mole numbers of all species in the multiphase -// object + void MultiPhase::getMoles(doublereal* molNum) const { /* @@ -556,10 +473,7 @@ void MultiPhase::getMoles(doublereal* molNum) const } } } -//==================================================================================================================== -/// Set the species moles to the values in array \a n. The state -/// of each phase object is also updated to have the specified -/// composition and the mixture temperature and pressure. + void MultiPhase::setMoles(const doublereal* n) { if (!m_init) { @@ -590,7 +504,7 @@ void MultiPhase::setMoles(const doublereal* n) loc += nsp; } } -//==================================================================================================================== + void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) { vector_fp tmpMoles(m_nsp, 0.0); @@ -601,7 +515,7 @@ void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) } setMoles(DATA_PTR(tmpMoles)); } -//==================================================================================================================== + void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) { if (!m_init) { @@ -611,7 +525,7 @@ void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) m_press = Pres; updatePhases(); } -//==================================================================================================================== + void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal* n) { @@ -619,7 +533,7 @@ void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres, m_press = Pres; setMoles(n); } -//==================================================================================================================== + void MultiPhase::getElemAbundances(doublereal* elemAbundances) const { size_t eGlobal; @@ -628,8 +542,7 @@ void MultiPhase::getElemAbundances(doublereal* elemAbundances) const elemAbundances[eGlobal] = m_elemAbundances[eGlobal]; } } -//==================================================================================================================== -// Internal routine to calculate the element abundance vector + void MultiPhase::calcElemAbundances() const { size_t loc = 0; @@ -653,8 +566,7 @@ void MultiPhase::calcElemAbundances() const loc += nspPhase; } } -//==================================================================================================================== -/// The total mixture volume [m^3]. + doublereal MultiPhase::volume() const { int i; @@ -665,7 +577,7 @@ doublereal MultiPhase::volume() const } return sum; } -//==================================================================================================================== + doublereal MultiPhase::equilibrate(int XY, doublereal err, int maxsteps, int maxiter, int loglevel) { @@ -1001,7 +913,7 @@ void importFromXML(string infile, string id) } } #endif -//==================================================================================================================== + void MultiPhase::setTemperature(const doublereal T) { if (!m_init) { @@ -1025,14 +937,11 @@ void MultiPhase::checkElementArraySize(size_t mm) const } } -//==================================================================================================================== -// Name of element \a m. std::string MultiPhase::elementName(size_t m) const { return m_enames[m]; } -//==================================================================================================================== -// Index of element with name \a name. + size_t MultiPhase::elementIndex(const std::string& name) const { for (size_t e = 0; e < m_nel; e++) { @@ -1057,8 +966,6 @@ void MultiPhase::checkSpeciesArraySize(size_t kk) const } } -//==================================================================================================================== -// Name of species with global index \a k. std::string MultiPhase::speciesName(const size_t k) const { return m_snames[k]; @@ -1068,18 +975,18 @@ doublereal MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const { return m_atoms(mGlob, kGlob); } -//==================================================================================================================== + void MultiPhase::getMoleFractions(doublereal* const x) const { std::copy(m_moleFractions.begin(), m_moleFractions.end(), x); } -//==================================================================================================================== + std::string MultiPhase::phaseName(const size_t iph) const { const ThermoPhase* tptr = m_phase[iph]; return tptr->id(); } -//==================================================================================================================== + int MultiPhase::phaseIndex(const std::string& pName) const { std::string tmp; @@ -1092,12 +999,12 @@ int MultiPhase::phaseIndex(const std::string& pName) const } return -1; } -//==================================================================================================================== + doublereal MultiPhase::phaseMoles(const size_t n) const { return m_moles[n]; } -//==================================================================================================================== + void MultiPhase::setPhaseMoles(const size_t n, const doublereal moles) { m_moles[n] = moles; @@ -1107,20 +1014,17 @@ size_t MultiPhase::speciesPhaseIndex(const size_t kGlob) const { return m_spphase[kGlob]; } -//==================================================================================================================== + doublereal MultiPhase::moleFraction(const size_t kGlob) const { return m_moleFractions[kGlob]; } -//==================================================================================================================== bool MultiPhase::tempOK(const size_t p) const { return m_temp_OK[p]; } -//==================================================================================================================== -/// Update the locally-stored species mole fractions. void MultiPhase::uploadMoleFractionsFromPhases() { size_t ip, loc = 0; @@ -1131,19 +1035,7 @@ void MultiPhase::uploadMoleFractionsFromPhases() } calcElemAbundances(); } -//==================================================================================================================== -//------------------------------------------------------------- -// -// protected methods -// -//------------------------------------------------------------- - - -/// synchronize the phase objects with the mixture state. This -/// method sets each phase to the mixture temperature and -/// pressure, and sets the phase mole fractions based on the -/// mixture mole numbers. void MultiPhase::updatePhases() const { size_t p, nsp, loc = 0; @@ -1159,6 +1051,4 @@ void MultiPhase::updatePhases() const } } } -//==================================================================================================================== } -