From 120dca66c7d86cb402326247962fc17eaf7b20f7 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Wed, 12 Nov 2003 18:58:15 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/python/Cantera/ctml_writer.py | 13 +- Cantera/python/examples/flame1.py | 3 +- Cantera/src/ArrayViewer.h | 4 - Cantera/src/BandMatrix.h | 4 - Cantera/src/ChemEquil.h | 6 - Cantera/src/Constituents.cpp | 1 - Cantera/src/Constituents.h | 194 ++--- Cantera/src/Group.h | 3 - Cantera/src/Phase.cpp | 9 + Cantera/src/Phase.h | 6 + Cantera/src/ReactionPath.cpp | 4 - Cantera/src/ReactionPath.h | 7 - Cantera/src/ReactionStoichMgr.cpp | 80 ++- Cantera/src/ReactionStoichMgr.h | 282 +++++--- Cantera/src/StoichManager.h | 975 ++++++++++++++------------ Cantera/src/config.h | 58 +- Cantera/src/ct_defs.h | 2 +- Cantera/src/misc.cpp | 46 +- Cantera/src/units.h | 1 + 19 files changed, 913 insertions(+), 785 deletions(-) diff --git a/Cantera/python/Cantera/ctml_writer.py b/Cantera/python/Cantera/ctml_writer.py index bc6d3538e..2abe48c7c 100644 --- a/Cantera/python/Cantera/ctml_writer.py +++ b/Cantera/python/Cantera/ctml_writer.py @@ -781,6 +781,9 @@ class phase(writer): self._sp = [] self._rx = [] self._options = options + self.debug = 0 + if 'debug' in options: + self.debug = 1 #-------------------------------- # process species @@ -943,12 +946,15 @@ class ideal_gas(phase): transport = 'None', initial_state = None, options = []): - + + print 'ig options = ',options phase.__init__(self, name, 3, elements, species, reactions, initial_state, options) self._pure = 0 self._kin = kinetics self._tr = transport + if self.debug: + print 'Read ideal_gas entry '+self._name def build(self, p): @@ -1198,7 +1204,10 @@ if __name__ == "__main__": # $Revision$ # $Date$ # $Log$ -# Revision 1.22 2003-11-01 04:48:20 dggoodwin +# Revision 1.23 2003-11-12 18:58:15 dggoodwin +# *** empty log message *** +# +# Revision 1.22 2003/11/01 04:48:20 dggoodwin # added capability to have species names with embedded commas # # Revision 1.21 2003/10/14 06:48:07 dggoodwin diff --git a/Cantera/python/examples/flame1.py b/Cantera/python/examples/flame1.py index df1895edf..c7fa58b64 100755 --- a/Cantera/python/examples/flame1.py +++ b/Cantera/python/examples/flame1.py @@ -16,6 +16,7 @@ tburner = 373.0 # burner temperature mdot = 0.06 # kg/m^2/s rxnmech = 'h2o2.cti' # reaction mechanism file +mix = 'ohmech' # gas mixture model comp = 'H2:1.8, O2:1, AR:7' # premixed gas composition # The solution domain is chosen to be 50 cm, and a point very near the @@ -40,7 +41,7 @@ refine_grid = 1 # 1 to enable refinement, 0 to # This object will be used to evaluate all thermodynamic, kinetic, # and transport properties # -gas = IdealGasMix(rxnmech) +gas = IdealGasMix(rxnmech, mix) # set its state to that of the unburned gas at the burner gas.setState_TPX(tburner, p, comp) diff --git a/Cantera/src/ArrayViewer.h b/Cantera/src/ArrayViewer.h index 906e5a990..d1cb8647d 100755 --- a/Cantera/src/ArrayViewer.h +++ b/Cantera/src/ArrayViewer.h @@ -15,10 +15,6 @@ #ifndef CT_ARRAYVIEWER_H #define CT_ARRAYVIEWER_H -#include -#include -using namespace std; - #include "ct_defs.h" #include "ctexceptions.h" #include "stringUtils.h" diff --git a/Cantera/src/BandMatrix.h b/Cantera/src/BandMatrix.h index 776282d74..6ab4dc845 100755 --- a/Cantera/src/BandMatrix.h +++ b/Cantera/src/BandMatrix.h @@ -14,10 +14,6 @@ #ifndef CT_BANDMATRIX_H #define CT_BANDMATRIX_H -#include -#include -using namespace std; - #include "ct_defs.h" #include "ctlapack.h" #include "utilities.h" diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index 2d5c5c272..35423d4e9 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -16,12 +16,6 @@ #define CT_CHEM_EQUIL_H -// STL includes -#include -#include -using namespace std; - - // Cantera includes #include "ct_defs.h" #include "vec_functions.h" diff --git a/Cantera/src/Constituents.cpp b/Cantera/src/Constituents.cpp index 8a923e409..492f24fd0 100755 --- a/Cantera/src/Constituents.cpp +++ b/Cantera/src/Constituents.cpp @@ -19,7 +19,6 @@ #include "Constituents.h" #include "Elements.h" -#include using namespace std; namespace Cantera { diff --git a/Cantera/src/Constituents.h b/Cantera/src/Constituents.h index 9cc8075d5..9a53f5737 100755 --- a/Cantera/src/Constituents.h +++ b/Cantera/src/Constituents.h @@ -1,21 +1,22 @@ -/** - * @file Constituents.h - * Header file for class Constituents - * - * $Author$ + +/// @file Constituents.h +/// Header file for class Constituents + + +/* $Author$ * $Date$ * $Revision$ * * $Log$ - * Revision 1.4 2003-09-03 18:15:50 hkmoffa + * Revision 1.5 2003-11-12 18:58:17 dggoodwin + * *** empty log message *** + * + * Revision 1.4 2003/09/03 18:15:50 hkmoffa * Added a vector get for the atoms in a species. * * Revision 1.3 2003/07/21 16:02:53 hkmoffa * Took out a double nested @name that gave a warning to doxygen * - * Revision 1.2 2003/06/27 14:19:16 dggoodwin - * *** empty log message *** - * * Revision 1.1.1.1 2003/04/14 17:57:51 dggoodwin * Initial import. * @@ -35,7 +36,6 @@ #include "ct_defs.h" using namespace std; -//#include "Elements.h" #include "SpeciesThermo.h" #include "ctexceptions.h" #include "stringUtils.h" @@ -45,24 +45,6 @@ namespace Cantera { class Elements; -#ifdef INCL_DEPRECATED_METHODS - - /** - * Structure returned by method species() - * @param name species name - * @param atoms vector of element atom numbers - * @param phase flag specifying phase - * @param charge electric charge - * @param molecularWeight molecular weight - */ - struct SpeciesData { - string name; - vector_fp atoms; - int phase; - doublereal charge; - doublereal molecularWeight; - }; -#endif /************** DEFINITIONS OF ERRORS *****************************/ @@ -73,86 +55,116 @@ namespace Cantera { " outside valid range of 0 to " + int2str(kmax-1)) {} }; + /******************************************************************/ + + + /// Class Constituents manages a set of elements and + /// species. Class Constituents is designed to provide information + /// about the elements and species in a phase - names, index + /// numbers (location in arrays), atomic or molecular weights, + /// etc. No computations are performed by the methods of this + /// class. The set of elements must include all those that compose + /// the species, but may include additional elements. The species + /// all must belong to the same phase. - /** - * Class Constituents manages a set of elements and species. The - * set of elements must include all those that compose the - * species, but may include additional elements. The species all - * must belong to the same phase. - */ class Constituents { public: - Constituents(Elements* ptr_Elements = 0); - ~Constituents(); + /// Constructor. + Constituents(Elements* ptr_Elements = 0); - /// Atomic weight of element m. - doublereal atomicWeight(int m) const; - /// vector of atomic weights - const array_fp& atomicWeights() const; - /// Number of elements. - int nElements() const; + /// Destructor. + ~Constituents(); -#ifdef INCL_DEPRECATED_METHODS - /** - * Returns an ElementData struct that contains the - * parameters for element number m. - */ - ElementData element(int m) const { - return m_Elements->element(m); - } -#endif - /** - * @name Adding Elements and Species - * These methods are used to add new elements or species. - * They are not usually called by user programs. - */ - void addElement(const string& symbol, doublereal weight); - void addElement(const XML_Node& e); - void addUniqueElement(const string& symbol, doublereal weight); - void addUniqueElement(const XML_Node& e); - /** - * Prohibit addition of more elements, and prepare to add - * species. - */ - void freezeElements(); - /// True if freezeElements has been called. - bool elementsFrozen(); - /** - * Index of element named 'name'. The index is an integer - * assigned to each element in the order it was added, - * beginning with 0 for the first element. If 'name' is not - * the name of an element in the set, then the value -1 is - * returned. - */ - int elementIndex(string name) const; - /** - * Name of the element with index m. @param m Element - * index. If m < 0 or m >= nElements() an exception is thrown. - */ - string elementName(int m) const; - /** - * Returns a read-only reference to the vector of element names. - */ - const vector& elementNames() const; + /// @name Element Information + //@{ + + /// Name of the element with index m. @param m Element + /// index. If m < 0 or m >= nElements() an exception is thrown. + string elementName(int m) const; - /** - * Returns the Number of species in the phase - */ + + /// Index of element named 'name'. The index is an integer + /// assigned to each element in the order it was added, + /// beginning with 0 for the first element. If 'name' is not + /// the name of an element in the set, then the value -1 is + /// returned. + int elementIndex(string name) const; + + + /// Atomic weight of element m. + doublereal atomicWeight(int m) const; + + + /// Return a read-only reference to the vector of element names. + const vector& elementNames() const; + + + /// Return a read-only reference to the vector of atomic weights. + const array_fp& atomicWeights() const; + + + /// Number of elements. + int nElements() const; + + //@} + + + + /// @name Adding Elements and Species + /// These methods are used to add new elements or species. + /// These are not usually called by user programs. + /// + /// Since species are checked to insure that they are only + /// composed of declared elements, it is necessary to first + /// add all elements before adding any species. + + //@{ + + /// Add an element. + /// @param symbol Atomic symbol string. + /// @param weight Atomic mass in amu. + void addElement(const string& symbol, doublereal weight); + + /// Add an element from an XML specification. + void addElement(const XML_Node& e); + + void addUniqueElement(const string& symbol, doublereal weight); + + void addUniqueElement(const XML_Node& e); + + + /// Prohibit addition of more elements, and prepare to add + /// species. + void freezeElements(); + + /// True if freezeElements has been called. + bool elementsFrozen(); + + //@} + + /// Returns the number of species in the phase int nSpecies() const { return m_kk; } + /// Molecular weight of species k. doublereal molecularWeight(int k) const; + + /// Molar mass. Preferred name for molecular weight. + doublereal molarMass(int k) const { + return molecularWeight(k); + } + /** * Return a const reference to the vector of molecular weights * of the species */ const array_fp& molecularWeights() const; - /** - * Electrical charge of one species k molecule, divided by - * \f$ e = 1.602 \times 10^{-19}\f$ Coulombs. - */ + + /// Electrical charge of one species k molecule, divided by + /// the magnitude of the electron charge ( \f$ e = 1.602 + /// \times 10^{-19}\f$ Coulombs). Dimensionless. doublereal charge(int k) const; /** diff --git a/Cantera/src/Group.h b/Cantera/src/Group.h index f094fe3c6..b62aa4cf4 100755 --- a/Cantera/src/Group.h +++ b/Cantera/src/Group.h @@ -12,9 +12,6 @@ #ifndef CT_RXNPATH_GROUP #define CT_RXNPATH_GROUP -#include -using namespace std; - #include "ct_defs.h" namespace Cantera { diff --git a/Cantera/src/Phase.cpp b/Cantera/src/Phase.cpp index 1c96075d0..32f8efda5 100755 --- a/Cantera/src/Phase.cpp +++ b/Cantera/src/Phase.cpp @@ -150,10 +150,19 @@ namespace Cantera { /** * Copy the vector of molecular weights into array weights. + * @deprecated */ void Phase::getMolecularWeights(int iwt, doublereal* weights) { const array_fp& mw = Constituents::molecularWeights(); copy(mw.begin(), mw.end(), weights); + } + + /** + * Copy the vector of molecular weights into array weights. + */ + void Phase::getMolecularWeights(doublereal* weights) { + const array_fp& mw = Constituents::molecularWeights(); + copy(mw.begin(), mw.end(), weights); } /** diff --git a/Cantera/src/Phase.h b/Cantera/src/Phase.h index 09e6cca05..ba0ae2b0d 100755 --- a/Cantera/src/Phase.h +++ b/Cantera/src/Phase.h @@ -31,6 +31,7 @@ namespace Cantera { * These classes are used to represent phases of matter. */ + /** * Base class for phases of matter. Class Phase derives from both * Constituents and State. In addition to the methods of those two @@ -130,6 +131,11 @@ namespace Cantera { */ void getMolecularWeights(int iwt, doublereal* weights); + /** + * Copy the vector of molecular weights into array weights. + */ + void getMolecularWeights(doublereal* weights); + /** * Return a const reference to the internal vector of * molecular weights. diff --git a/Cantera/src/ReactionPath.cpp b/Cantera/src/ReactionPath.cpp index cd264d5ad..a7cb79498 100755 --- a/Cantera/src/ReactionPath.cpp +++ b/Cantera/src/ReactionPath.cpp @@ -16,10 +16,6 @@ #pragma warning(disable:4503) #endif -// STL includes -#include -#include - #include "ReactionPath.h" #include "Kinetics.h" #include "reaction_defs.h" diff --git a/Cantera/src/ReactionPath.h b/Cantera/src/ReactionPath.h index 5e51dfca2..4de70cb83 100755 --- a/Cantera/src/ReactionPath.h +++ b/Cantera/src/ReactionPath.h @@ -14,13 +14,6 @@ #ifndef CT_RXNPATH_H #define CT_RXNPATH_H -// STL includes -#include -#include -#include -#include -using namespace std; - // Cantera includes #include "ct_defs.h" #include "DenseMatrix.h" diff --git a/Cantera/src/ReactionStoichMgr.cpp b/Cantera/src/ReactionStoichMgr.cpp index a776d4990..37adee654 100644 --- a/Cantera/src/ReactionStoichMgr.cpp +++ b/Cantera/src/ReactionStoichMgr.cpp @@ -10,43 +10,57 @@ namespace Cantera { - ReactionStoichMgr:: - ReactionStoichMgr() { - m_reactants = new StoichManagerN; - m_revproducts = new StoichManagerN; - m_irrevproducts = new StoichManagerN; - } + // create stoichiometry managers for the reactants of all reactions, + // for the products of the reversible reactions, and for the + // products of the irreversible reactions. + ReactionStoichMgr:: + ReactionStoichMgr() { + m_reactants = new StoichManagerN; + m_revproducts = new StoichManagerN; + m_irrevproducts = new StoichManagerN; + } - ReactionStoichMgr::~ReactionStoichMgr() { - delete m_reactants; - delete m_revproducts; - delete m_irrevproducts; - } + // delete the three stoichiometry managers + ReactionStoichMgr::~ReactionStoichMgr() { + delete m_reactants; + delete m_revproducts; + delete m_irrevproducts; + } - void ReactionStoichMgr:: - add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible) { - vector_fp forder(reactants.size(), 1.0); - add(rxn, reactants, products, reversible, forder); - } - void ReactionStoichMgr:: - add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible, const vector_fp& fwdOrder) { - m_reactants->add(rxn, reactants, fwdOrder); - if (reversible) - m_revproducts->add(rxn, products); - else - m_irrevproducts->add(rxn, products); - } + void ReactionStoichMgr:: + add(int rxn, const vector_int& reactants, const vector_int& products, + bool reversible) { + vector_fp forder(reactants.size(), 1.0); + add(rxn, reactants, products, reversible, forder); + } - void ReactionStoichMgr:: - getCreationRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c) { - fill(c, c + nsp, 0.0); - m_revproducts->incrementSpecies(ropf, c); - m_irrevproducts->incrementSpecies(ropf, c); - m_reactants->incrementSpecies(ropr, c); - } + + void ReactionStoichMgr:: + add(int rxn, const vector_int& reactants, const vector_int& products, + bool reversible, const vector_fp& fwdOrder) { + + // add the reactants with the specified forward order + m_reactants->add(rxn, reactants, fwdOrder); + + // depending on whether the reversible flag is set or not, add the + // products either to the reversible or irreversible product + // stoichiometry manager. + if (reversible) + m_revproducts->add(rxn, products); + else + m_irrevproducts->add(rxn, products); + } + + + void ReactionStoichMgr:: + getCreationRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c) { + // zero out the target array + fill(c, c + nsp, 0.0); + m_revproducts->incrementSpecies(ropf, c); + m_irrevproducts->incrementSpecies(ropf, c); + m_reactants->incrementSpecies(ropr, c); + } void ReactionStoichMgr:: getDestructionRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* d) { diff --git a/Cantera/src/ReactionStoichMgr.h b/Cantera/src/ReactionStoichMgr.h index a2bece857..9b9e23a0c 100644 --- a/Cantera/src/ReactionStoichMgr.h +++ b/Cantera/src/ReactionStoichMgr.h @@ -1,6 +1,10 @@ /** * @file ReactionStoichMgr.h * + * Header file declaring class ReactionStoichMgr. + */ + +/* * $Author$ * $Revision$ * $Date$ @@ -9,113 +13,217 @@ #ifndef CT_RXN_STOICH #define CT_RXN_STOICH - #include "ct_defs.h" namespace Cantera { - class StoichManagerN; + class StoichManagerN; + + /** + * Reaction mechanism stoichiometry manager. This is an internal class used + * by kinetics manager classes, and is not meant for direct use in + * user programs. + * + * Class ReactionStoichMgr handles the calculation of quantities involving + * the stoichiometry of a set of reactions. The reactions must have integer + * stoichiometric coefficients. Specifically, its methods compute + * - species creation rates + * - species destruction rates + * - species net production rates + * - the change in molar species properties in the reactions + * - concentration products + * + * To use this class, method 'add' is first used to add each reaction. + * Once all reactions have been added, the methods that compute various + * quantities may be called. + * + * The nomenclature used below to document the methods is as follows. + * - \f$ N_r \f$ + * - Integer reactant stoichiometric coefficient matrix. The (k,i) + * element of this matrix is the stoichiometric coefficient of + * species \i k as a reactant in reaction \i i. + * - \f$ N_p \f$ + * - Integer product stoichiometric coefficient matrix. The (k,i) + * element of this matrix is the stoichiometric coefficient of + * species \i k as a product in reaction \i i. + * - \f$ Q_{\rm fwd} \f$ + * - Vector of length I of forward rates of progress. + * - \f$ Q_{\rm rev} \f$ + * - Vector of length I of reverse rates of progress. + * - \f$ C \f$ + * - Vector of K species creation rates. + * - \f$ D \f$ + * - Vector of K species destruction rates. + * - \f$ W = C - D \f$ + * - Vector of K species net production rates. + * + */ + class ReactionStoichMgr { + + public: + + /// Constructor. + ReactionStoichMgr(); + + /// Destructor. + virtual ~ReactionStoichMgr(); /** - * This class handles calculations involving reaction stoichiometry. + * Add a reaction with mass-action kinetics. Vectors + * 'reactants' and 'products' contain the integer species + * indices of the reactants and products, respectively. Note + * that if more than one molecule of a given species is + * involved in the reaction, then its index is repeated. + * + * For example, suppose a reaction mechanism involves the + * species N2, O2, O, N, NO. N2 is assigned index number 0, O2 + * number 1, and so on through NO with number 4. Then the + * representation of the following reactions is as shown here. + * + * - N + O = NO + * - reactants: (3, 2) + * - products: (4) + * + * - O + O = O2 + * - reactants: (2, 2) [ note repeated index ] + * - products: (1) + * + * @param rxn Reaction number. This number will be used as the index into the + * rate of progess vector in the methods below. + * @param reactants vector of integer reactant indices + * @param products vector of integer product indices + * @param reversible true if the reaction is reversible, false otherwise */ - class ReactionStoichMgr { + void add(int rxn, const vector_int& reactants, const vector_int& products, + bool reversible); - public: + /** + * Add a reaction with specified, possibly non-integral, reaction orders. + * @param rxn Reaction number + * @param reactants vector of integer reactant indices + * @param products vector of integer product indices + * @param reversible true if the reaction is reversible, false otherwise. + * If the reaction is reversible, its reverse rate will be computed from + * the reaction stoichiometry. + * @param fwdOrder reaction orders for the reactants. This vector must + * be the same length as 'reactants,' and the reaction orders are for the + * species with index in the corresponding location in 'reactants.' + * + */ + void add(int rxn, const vector_int& reactants, const vector_int& products, + bool reversible, const vector_fp& fwdOrder); - ReactionStoichMgr(); - virtual ~ReactionStoichMgr(); - /** - * Add a reaction with specified, possibly non-integral, reaction orders. - * @param rxn Reaction number - * @param reactants vector of integer reactant indices - * @param reactants vector of integer product indices - * @param reversible true if the reaction is reversible, false otherwise - * @param fwdOrder reaction orders for the reactants. This vector must - * be the same length as 'reactants,' and the reaction orders are for the - * species with index in the corresponding location in 'reactants.' - */ - void add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible, const vector_fp& fwdOrder); - /** - * Add a reaction with mass-action kinetics. - * @param rxn Reaction number - * @param reactants vector of integer reactant indices - * @param reactants vector of integer product indices - * @param reversible true if the reaction is reversible, false otherwise - */ - void add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible); - /** - * Given the arrays of the forward and reverse rates of progress for all reactions, - * compute the species creation rates and return them in array c. - */ - void getCreationRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c); + /** + * Species creation rates. + * Given the arrays of the forward and reverse rates of + * progress for all reactions, compute the species creation + * rates, given by + * \f[ + * C = N_p Q_f + N_r Q_r. + * \f] + */ + void getCreationRates(int nSpecies, + const doublereal* fwdRatesOfProgress, + const doublereal* revRatesOfProgress, + doublereal* creationRates); - /** - * Given the arrays of the forward and reverse rates of progress for all reactions, - * compute the species destruction rates and return them in array d. - */ - void getDestructionRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* d); - /** - * Given the array of the net rates of progress for all reactions, - * compute the species net production rates and return them in array w. - */ - void getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w); + /** + * Species destruction rates. + * Given the arrays of the forward and reverse rates of + * progress for all reactions, compute the species destruction + * rates, given by + * \f[ + * D = N_r Q_f + N_p Q_r, + * \f] + * Note that the stoichiometric coefficient matrices are very sparse, integer + * matrices. + */ + void getDestructionRates(int nSpecies, + const doublereal* fwdRatesOfProgress, + const doublereal* revRatesOfProgress, + doublereal* destructionRates); - /** - * Given an array of species properties 'g', return in array 'dg' the change in this quantity - * in the reactions. Array 'g' must have a length at least as great - * as the number of species, and array 'dg' must have a length - * as great as the total number of reactions. - */ - void getReactionDelta(int nr, const doublereal* g, doublereal* dg); - /** - * Given an array of species properties 'g', return in array - * 'dg' the change in this quantity in the reversible - * reactions. Array 'g' must have a length at least as great - * as the number of species, and array 'dg' must have a length - * as great as the total number of reactions. This method - * only computes 'dg' for the reversible reactions, and the - * entries of 'dg' for the irreversible reactions are - * unaltered. This is primarily designed for use in - * calculating reveerse rate coefficients from thermochemistry - * for reversible reactions. - */ - void getRevReactionDelta(int nr, const doublereal* g, doublereal* dg); + /** + * Given the array of the net rates of progress for all + * reactions, compute the species net production rates and + * return them in array w. + */ + /** + * Species net production rates. + * Given the array of the net rates of + * progress for all reactions, compute the species net production + * rates, given by + * \f[ + * W = (N_r - N_p) Q_{\rm net}, + * \f] + */ + void getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w); - /** - * Given an array of concentrations C, multiply the entries in array R by - * the concentration products for the reactants: - * \f[ - * R_i = R_i * \prod_k C_k^{o_{k,i}} - * \f] - * Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i. - */ - void multiplyReactants(const doublereal* C, doublereal* R); - /** - * Given an array of concentrations C, multiply the entries in array R by - * the concentration products for the products: - * \f[ - * R_i = R_i * \prod_k C_k^{\nu^{(p)}_{k,i}} - * \f] - * Here \f$ \nu^{(p)}_{k,i} \f$ is the product stoichiometric coefficient - * of species k in reaction i. - */ - void multiplyRevProducts(const doublereal* c, doublereal* r); - protected: + /** + * Change of a molar species property in a reaction. Given an + * array of species properties 'g', return in array 'dg' the + * change in this quantity in the reactions. Array 'g' must + * have a length at least as great as the number of species, + * and array 'dg' must have a length as great as the total + * number of reactions. + */ + void getReactionDelta(int nReactions, + const doublereal* g, + doublereal* dg); - StoichManagerN* m_reactants; - StoichManagerN* m_revproducts; - StoichManagerN* m_irrevproducts; - }; + + /** + * Given an array of species properties 'g', return in array + * 'dg' the change in this quantity in the reversible + * reactions. Array 'g' must have a length at least as great + * as the number of species, and array 'dg' must have a length + * as great as the total number of reactions. This method + * only computes 'dg' for the reversible reactions, and the + * entries of 'dg' for the irreversible reactions are + * unaltered. This is primarily designed for use in + * calculating reveerse rate coefficients from thermochemistry + * for reversible reactions. + */ + void getRevReactionDelta(int nr, const doublereal* g, doublereal* dg); + + + /** + * Given an array of concentrations C, multiply the entries in array R by + * the concentration products for the reactants: + * \f[ + * R_i = R_i * \prod_k C_k^{o_{k,i}} + * \f] + * Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i. + */ + void multiplyReactants(const doublereal* C, doublereal* R); + + + /** + * Given an array of concentrations C, multiply the entries in array R by + * the concentration products for the products: + * \f[ + * R_i = R_i * \prod_k C_k^{\nu^{(p)}_{k,i}} + * \f] + * Here \f$ \nu^{(p)}_{k,i} \f$ is the product stoichiometric coefficient + * of species k in reaction i. + */ + void multiplyRevProducts(const doublereal* c, doublereal* r); + + + protected: + + StoichManagerN* m_reactants; + StoichManagerN* m_revproducts; + StoichManagerN* m_irrevproducts; + + }; } #endif diff --git a/Cantera/src/StoichManager.h b/Cantera/src/StoichManager.h index 7461a761f..440215e4c 100755 --- a/Cantera/src/StoichManager.h +++ b/Cantera/src/StoichManager.h @@ -20,507 +20,560 @@ using namespace std; namespace Cantera { - /** - * @defgroup Stoichiometry Stoichiometry - * - * Note: these classes are designed for internal use in class ReactionStoichManager. - * - * Operations on reactions that require knowing the reaction - * stoichiometry. This module consists of class StoichManager, and - * classes C1, C2, and C3. Classes C1, C2, and C3 handle operations - * involving one, two, or three species, respectively, in a - * reaction. Instances are instantiated with a reaction number, and n - * species numbers (n = 1 for C1, etc.). All three classes have the - * same interface. - * - * These classes are designed for use by StoichManager, and the - * operations implemented are those needed to efficiently compute - * quantities such as rates of progress, species production rates, - * reaction thermochemistry, etc. The compiler will inline these - * methods into the body of the corresponding StoichManager method, - * and so there is no performance penalty (unless inlining is turned - * off). - * - * To describe the methods, consider class C3 and suppose an instance - * is created with reaction number irxn and species numbers k0, k1, - * and k2. - * - * - multiply(in, out) : out[irxn] is multiplied by - * in[k0] * in[k1] * in[k2] - * - * - power(in, out) : out[irxn] is multiplied by - * (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2) - * - * - incrementReaction(in, out) : out[irxn] is incremented by - * in[k0] + in[k1] + in[k2] - * - * - decrementReaction(in, out) : out[irxn] is decremented by - * in[k0] + in[k1] + in[k2] - * - * - incrementSpecies(in, out) : out[k0], out[k1], and out[k2] - * are all incremented by in[irxn] - * - * - decrementSpecies(in, out) : out[k0], out[k1], and out[k2] - * are all decremented by in[irxn] - * - * The function multiply() is usually used when evaluating the - * forward and reverse rates of progress of reactions. - * The rate constants are usually loaded into out[]. Then - * multply() is called to add in the dependence of the - * species concentrations to yield a forward and reverse rop. - * - * The function incrementSpecies() and its cousin decrementSpecies() - * is used to translate from rates of progress to species production - * rates. The vector in[] is preloaed with the rates of progess of - * all reactions. Then incrementSpecies() is called to - * increment the species production vector, out[], with the rates - * of progress. - * - * The functions incrementReaction() and decrementReaction() are - * used to find the standard state equilibrium constant for - * a reaction. Here, output[] is a vector of length - * number of reactions, usually the standard gibbs free energies - * of reaction, while input, usually the standard state - * gibbs free energies of species, is a vector of length number of - * species. - * - * Note the stoichiometric coefficient for a species in a reaction - * is handled by always assuming it is equal to one and then - * treating reactants and products for a reaction separately. - * Bimolecular reactions involving the identical species are - * treated as involving separate species. - * - * @internal This class should be upgraded to include cases where - * real stoichiometric coefficients are used. Shouldn't be that - * hard to do, and they occur in engineering simulations with some - * regularity. - * - */ + /** + * @defgroup Stoichiometry Stoichiometry + * + * Note: these classes are designed for internal use in class + * ReactionStoichManager. + * + * Operations on reactions that require knowing the reaction + * stoichiometry. This module consists of class StoichManager, and + * classes C1, C2, and C3. Classes C1, C2, and C3 handle operations + * involving one, two, or three species, respectively, in a + * reaction. Instances are instantiated with a reaction number, and n + * species numbers (n = 1 for C1, etc.). All three classes have the + * same interface. + * + * These classes are designed for use by StoichManager, and the + * operations implemented are those needed to efficiently compute + * quantities such as rates of progress, species production rates, + * reaction thermochemistry, etc. The compiler will inline these + * methods into the body of the corresponding StoichManager method, + * and so there is no performance penalty (unless inlining is turned + * off). + * + * To describe the methods, consider class C3 and suppose an instance + * is created with reaction number irxn and species numbers k0, k1, + * and k2. + * + * - multiply(in, out) : out[irxn] is multiplied by + * in[k0] * in[k1] * in[k2] + * + * - power(in, out) : out[irxn] is multiplied by + * (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2) + * + * - incrementReaction(in, out) : out[irxn] is incremented by + * in[k0] + in[k1] + in[k2] + * + * - decrementReaction(in, out) : out[irxn] is decremented by + * in[k0] + in[k1] + in[k2] + * + * - incrementSpecies(in, out) : out[k0], out[k1], and out[k2] + * are all incremented by in[irxn] + * + * - decrementSpecies(in, out) : out[k0], out[k1], and out[k2] + * are all decremented by in[irxn] + * + * The function multiply() is usually used when evaluating the + * forward and reverse rates of progress of reactions. + * The rate constants are usually loaded into out[]. Then + * multply() is called to add in the dependence of the + * species concentrations to yield a forward and reverse rop. + * + * The function incrementSpecies() and its cousin decrementSpecies() + * is used to translate from rates of progress to species production + * rates. The vector in[] is preloaed with the rates of progess of + * all reactions. Then incrementSpecies() is called to + * increment the species production vector, out[], with the rates + * of progress. + * + * The functions incrementReaction() and decrementReaction() are + * used to find the standard state equilibrium constant for + * a reaction. Here, output[] is a vector of length + * number of reactions, usually the standard gibbs free energies + * of reaction, while input, usually the standard state + * gibbs free energies of species, is a vector of length number of + * species. + * + * Note the stoichiometric coefficient for a species in a reaction + * is handled by always assuming it is equal to one and then + * treating reactants and products for a reaction separately. + * Bimolecular reactions involving the identical species are + * treated as involving separate species. + * + * @internal This class should be upgraded to include cases where + * real stoichiometric coefficients are used. Shouldn't be that + * hard to do, and they occur in engineering simulations with some + * regularity. + * + */ - /** - * Handles one molecule in a reaction. - * @ingroup Stoichiometry - * @internal - */ - class C1 { + /** + * Handles one species in a reaction. + * @ingroup Stoichiometry + * @internal + */ + class C1 { - public: + public: - C1( int rxn = 0, int ic0 = 0, doublereal order = 1.0 ) - : m_rxn (rxn), m_ic0 (ic0), m_order(order) {} + C1( int rxn = 0, int ic0 = 0, doublereal order = 1.0 ) + : m_rxn (rxn), m_ic0 (ic0), m_order(order) {} - int data(vector& ic) { - ic.resize(3); - ic[0] = m_ic0; - return m_rxn; - } + int data(vector& ic) { + ic.resize(3); + ic[0] = m_ic0; + return m_rxn; + } - void multiply(const doublereal* input, doublereal* output) const { - *(output + m_rxn) *= *(input + m_ic0); - } - void power(const doublereal* input, doublereal* output) const { - output[m_rxn] *= pow(input[m_ic0], m_order); - } - void incrementSpecies(const doublereal* input, - doublereal* output) const { - *(output + m_ic0) += *(input + m_rxn); - } - void decrementSpecies(const doublereal* input, - doublereal* output) const { - *(output + m_ic0) -= *(input + m_rxn); - } - void incrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) += *(input + m_ic0); - } - void decrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) -= *(input + m_ic0); - } - private: - int m_rxn, m_ic0; - doublereal m_order; - }; + void multiply(const doublereal* input, doublereal* output) const { + *(output + m_rxn) *= *(input + m_ic0); + } + void power(const doublereal* input, doublereal* output) const { + output[m_rxn] *= pow(input[m_ic0], m_order); + } + void incrementSpecies(const doublereal* input, + doublereal* output) const { + *(output + m_ic0) += *(input + m_rxn); + } + void decrementSpecies(const doublereal* input, + doublereal* output) const { + *(output + m_ic0) -= *(input + m_rxn); + } + void incrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) += *(input + m_ic0); + } + void decrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) -= *(input + m_ic0); + } + private: + int m_rxn, m_ic0; + doublereal m_order; + }; + /** + * Handles two species in a single reaction. + * @ingroup Stoichiometry + */ + class C2 { + public: + C2( int rxn = 0, int ic0 = 0, int ic1 = 0, + doublereal order0 = 1.0, doublereal order1 = 1.0 ) + : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), + m_order0(order0), m_order1(order1) {} + + int data(vector& ic) { + ic.resize(2); + ic[0] = m_ic0; + ic[1] = m_ic1; + return m_rxn; + } + + void multiply(const doublereal* input, doublereal* output) const { + output[m_rxn] *= input[m_ic0] * input[m_ic1]; + } + + void power(const doublereal* input, doublereal* output) const { + output[m_rxn] *= pow(input[m_ic0],m_order0) * + pow(input[m_ic1],m_order1); + } + + void incrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + output[m_ic0] += x; + output[m_ic1] += x; + } + + void decrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + output[m_ic0] -= x; + output[m_ic1] -= x; + } + + void incrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1); + } + + void decrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1)); + } + + private: + /** - * Handles two species in a single reaction. - * @ingroup Stoichiometry + * Reaction index -> index into the ROP vector */ - class C2 { - public: - C2( int rxn = 0, int ic0 = 0, int ic1 = 0, - doublereal order0 = 1.0, doublereal order1 = 1.0 ) - : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), - m_order0(order0), m_order1(order1) {} + int m_rxn; - int data(vector& ic) { - ic.resize(2); - ic[0] = m_ic0; - ic[1] = m_ic1; - return m_rxn; - } - - void multiply(const doublereal* input, doublereal* output) const { - output[m_rxn] *= input[m_ic0] * input[m_ic1]; - } - void power(const doublereal* input, doublereal* output) const { - output[m_rxn] *= pow(input[m_ic0],m_order0) * - pow(input[m_ic1],m_order1); - } - void incrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - output[m_ic0] += x; - output[m_ic1] += x; - } - void decrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - output[m_ic0] -= x; - output[m_ic1] -= x; - } - void incrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1); - } - void decrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1)); - } - private: - /** - * Reaction index -> index into the ROP vector - */ - int m_rxn; - /** - * Species indecise -> index into the species vector for the two species. - */ - int m_ic0, m_ic1; - doublereal m_order0, m_order1; - }; + /** + * Species indecise -> index into the species vector for the + * two species. + */ + int m_ic0, m_ic1; + doublereal m_order0, m_order1; + }; - /** - * Handles three species in a reaction. - * @ingroup Stoichiometry - */ - class C3 { - public: - C3( int rxn = 0, int ic0 = 0, int ic1 = 0, int ic2 = 0, - doublereal order0 = 1.0, doublereal order1 = 1.0, - doublereal order2 = 1.0) - : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), m_ic2 (ic2), - m_order0(order0), m_order1(order1), m_order2(order2) {} + /** + * Handles three species in a reaction. + * @ingroup Stoichiometry + */ + class C3 { + public: + C3( int rxn = 0, int ic0 = 0, int ic1 = 0, int ic2 = 0, + doublereal order0 = 1.0, doublereal order1 = 1.0, + doublereal order2 = 1.0) + : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), m_ic2 (ic2), + m_order0(order0), m_order1(order1), m_order2(order2) {} - int data(vector& ic) { - ic.resize(3); - ic[0] = m_ic0; - ic[1] = m_ic1; - ic[2] = m_ic2; - return m_rxn; - } + int data(vector& ic) { + ic.resize(3); + ic[0] = m_ic0; + ic[1] = m_ic1; + ic[2] = m_ic2; + return m_rxn; + } - void multiply(const doublereal* input, doublereal* output) const { - *(output + m_rxn) *= (*(input + m_ic0)) * (*(input + m_ic1)) - * (*(input + m_ic2)); - } - void power(const doublereal* input, doublereal* output) const { - output[m_rxn] *= pow(input[m_ic0],m_order0) * - pow(input[m_ic1],m_order1) * - pow(input[m_ic2],m_order2); - } - void incrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = *(input + m_rxn); - *(output + m_ic0) += x; - *(output + m_ic1) += x; - *(output + m_ic2) += x; - } - void decrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = *(input + m_rxn); - *(output + m_ic0) -= x; - *(output + m_ic1) -= x; - *(output + m_ic2) -= x; - } - void incrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1) - + *(input + m_ic2); - } - void decrementReaction(const doublereal* input, - doublereal* output) const { - *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1) - + *(input + m_ic2)); - } - private: - int m_rxn, m_ic0, m_ic1, m_ic2; - doublereal m_order0, m_order1, m_order2; - }; - - /** - * Handles any number of species in a reaction. - * @ingroup Stoichiometry - */ - class C_AnyN { - public: - C_AnyN() : m_rxn (-1) {} - - C_AnyN( int rxn, const vector_int& ic, const vector_fp& order) - : m_rxn (rxn) { - m_n = ic.size(); - m_ic.resize(m_n); - m_order.resize(m_n); - for (int n = 0; n < m_n; n++) { - m_ic[n] = ic[n]; - m_order[n] = order[n]; - } - } - - int data(vector& ic) { - ic.resize(m_n); - int n; - for (n = 0; n < m_n; n++) ic[n] = m_ic[n]; - return m_rxn; - } - - void power(const doublereal* input, doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] - *= pow(input[m_ic[n]],m_order[n]); - } - - void multiply(const doublereal* input, doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] *= input[m_ic[n]]; - } - - void incrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - for (int n = 0; n < m_n; n++) output[m_ic[n]] += x; - } - void decrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - for (int n = 0; n < m_n; n++) output[m_ic[n]] -= x; - } - void incrementReaction(const doublereal* input, - doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] += input[m_ic[n]]; - } - void decrementReaction(const doublereal* input, - doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] -= input[m_ic[n]]; - } - private: - int m_n, m_rxn, m_ic0, m_ic1, m_ic2; - vector_int m_ic; - vector_fp m_order; - }; + void multiply(const doublereal* input, doublereal* output) const { + *(output + m_rxn) *= (*(input + m_ic0)) * (*(input + m_ic1)) + * (*(input + m_ic2)); + } + void power(const doublereal* input, doublereal* output) const { + output[m_rxn] *= pow(input[m_ic0],m_order0) * + pow(input[m_ic1],m_order1) * + pow(input[m_ic2],m_order2); + } + void incrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = *(input + m_rxn); + *(output + m_ic0) += x; + *(output + m_ic1) += x; + *(output + m_ic2) += x; + } + void decrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = *(input + m_rxn); + *(output + m_ic0) -= x; + *(output + m_ic1) -= x; + *(output + m_ic2) -= x; + } + void incrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1) + + *(input + m_ic2); + } + void decrementReaction(const doublereal* input, + doublereal* output) const { + *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1) + + *(input + m_ic2)); + } + private: + int m_rxn, m_ic0, m_ic1, m_ic2; + doublereal m_order0, m_order1, m_order2; + }; - template + /** + * Handles any number of species in a reaction. + * @ingroup Stoichiometry + */ + class C_AnyN { + public: + C_AnyN() : m_rxn (-1) {} + + C_AnyN( int rxn, const vector_int& ic, const vector_fp& order) + : m_rxn (rxn) { + m_n = ic.size(); + m_ic.resize(m_n); + m_order.resize(m_n); + for (int n = 0; n < m_n; n++) { + m_ic[n] = ic[n]; + m_order[n] = order[n]; + } + } + + int data(vector& ic) { + ic.resize(m_n); + int n; + for (n = 0; n < m_n; n++) ic[n] = m_ic[n]; + return m_rxn; + } + + void power(const doublereal* input, doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] + *= pow(input[m_ic[n]],m_order[n]); + } + + void multiply(const doublereal* input, doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] *= input[m_ic[n]]; + } + + void incrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + for (int n = 0; n < m_n; n++) output[m_ic[n]] += x; + } + void decrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + for (int n = 0; n < m_n; n++) output[m_ic[n]] -= x; + } + void incrementReaction(const doublereal* input, + doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] += input[m_ic[n]]; + } + void decrementReaction(const doublereal* input, + doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] -= input[m_ic[n]]; + } + private: + int m_n, m_rxn, m_ic0, m_ic1, m_ic2; + vector_int m_ic; + vector_fp m_order; + }; + + + template inline static void _multiply(_InputIter __begin, _InputIter __end, - const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->multiply(input, output); - } - template + const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->multiply(input, output); + } + template inline static void _power(_InputIter __begin, _InputIter __end, - const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->power(input, output); - } - template + const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->power(input, output); + } + template inline static void _incrementSpecies(_InputIter __begin, - _InputIter __end, const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->incrementSpecies(input, output); - } - template + _InputIter __end, const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->incrementSpecies(input, output); + } + template inline static void _decrementSpecies(_InputIter __begin, - _InputIter __end, const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->decrementSpecies(input, output); - } - template + _InputIter __end, const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->decrementSpecies(input, output); + } + template inline static void _incrementReactions(_InputIter __begin, - _InputIter __end, const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->incrementReaction(input, output); - } - template + _InputIter __end, const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->incrementReaction(input, output); + } + template inline static void _decrementReactions(_InputIter __begin, - _InputIter __end, const Vec1& input, Vec2& output) { - for (; __begin != __end; ++__begin) - __begin->decrementReaction(input, output); + _InputIter __end, const Vec1& input, Vec2& output) { + for (; __begin != __end; ++__begin) + __begin->decrementReaction(input, output); + } + + + /** + * This class handles operations involving the stoichiometric + * coefficients on one side of a reaction (reactant or product) for + * a set of reactions comprising a reaction mechanism. This class is + * used by class ReactionStoichMgr, which contains three instances + * of this class (one to handle operations on the reactions, one for + * the products of reversible reactions, and one for the products of + * irreversible reactions). + * + * This class is designed for use with elementary reactions, or at + * least ones with integral stoichiometric coefficients. Let \f$ M(i) \f$ + * be the number of molecules on the product or reactant side of + * reaction number i. + * \f[ + * r_i = \sum_m^{M_i} s_{k_{m,i}} + * \f] + * To understand the operations performed by this class, let + * $N_{k,i}$ denote the stoichiometric coefficient of species k on + * one side (reactant or product) in reaction i. Then \b N is a sparse + * K by I matrix of stoichiometric coefficients. + * + * The following matrix operations may be carried out with a vector + * S of length K, and a vector R of length I: + * + * - \f$ S = S + N R\f$ (incrementSpecies) + * - \f$ S = S - N R\f$ (decrementSpecies) + * - \f$ R = R + N^T S \f$ (incrementReaction) + * - \f$ R = R - N^T S \f$ (deccrementReaction) + * + * The actual implementation, however, does not compute these + * quantities by matrix multiplication. A faster algorithm is used + * that makes use of the fact that the \b integer-valued N matrix is + * very sparse, and the non-zero terms are small positive integers. + * \f[ + * S_k = R_{i1} + \dots + R_{iM} + * \f] + * where M is the number of molecules, and $\f i(m) \f$ is the + * @ingroup Stoichiometry + */ + class StoichManagerN { + public: + + /** + * Constructor for the StoichManagerN class. + * + * @internal Consider adding defaulted entries here that supply + * the total number of reactions in the mechanism and the total + * number of species in the species list. Then, we could use those + * numbers to provide error checks during the construction of the + * object. Those numbers would also provide some clarity to the + * purpose and utility of this class. + * + * DGG - the problem is that the number of reactions and species + * are not known initially. + */ + StoichManagerN() {} + + /** + * Add a single reaction to the list of reactions that this + * stoichiometric manager object handles. + * + * This function is the same as the add() function below. However, + * the order of each species in the power list expression is + * set to one automatically. + */ + void add(int rxn, const vector_int& k) { + vector_fp order(k.size(), 1.0); + add(rxn, k, order); } - - class StoichManagerN { - public: - - /** - * Constructor for the StoichManagerN class. - * - * @internal Consider adding defaulted entries here that supply - * the total number of reactions in the mechanism and the total - * number of species in the species list. Then, we could use those - * numbers to provide error checks during the construction of the - * object. Those numbers would also provide some clarity to the - * purpose and utility of this class. - */ - StoichManagerN() {} - - /** - * Add a single reaction to the list of reactions that this - * stoichiometric manager object handles. - * - * This function is the same as the add() function below. However, - * the order of each species in the power list expression is - * set to one automatically. - */ - void add(int rxn, const vector_int& k) { - vector_fp order(k.size(), 1.0); - add(rxn, k, order); - } - - /** - * Add a single reaction to the list of reactions that this - * stoichiometric manager object handles. - * - * @param rxn Reaction index of the current reaction. This is used - * as an index into vectors which have length n_total_rxn. - * @param k This is a vector of integer values specifying the - * species indecises. The length of this vector species - * the number of different species in the description. - * The value of the entries are the species indices. - * These are used as indexes into vectors which have - * length n_total_species. - * @param order This is a vector of the same length as vector k. - * The order is used for the routine order(), which produces - * a power law expression involving the species vector. - */ - void add(int rxn, const vector_int& k, const vector_fp& order) { - m_n[rxn] = k.size(); - switch (k.size()) { - case 1: - m_loc[rxn] = m_c1_list.size(); - m_c1_list.push_back(C1(rxn, k[0], order[0])); - break; - case 2: - m_loc[rxn] = m_c2_list.size(); - m_c2_list.push_back(C2(rxn, k[0], k[1], order[0], order[1])); - break; - case 3: - m_loc[rxn] = m_c3_list.size(); - m_c3_list.push_back(C3(rxn, k[0], k[1], k[2], - order[0], order[1], order[2])); - break; - default: - m_loc[rxn] = m_cn_list.size(); - m_cn_list.push_back(C_AnyN(rxn, k, order)); - } - } + /** + * Add a single reaction to the list of reactions that this + * stoichiometric manager object handles. + * + * @param rxn Reaction index of the current reaction. This is used + * as an index into vectors which have length n_total_rxn. + * @param k This is a vector of integer values specifying the + * species indecises. The length of this vector species + * the number of different species in the description. + * The value of the entries are the species indices. + * These are used as indexes into vectors which have + * length n_total_species. + * @param order This is a vector of the same length as vector k. + * The order is used for the routine order(), which produces + * a power law expression involving the species vector. + */ + void add(int rxn, const vector_int& k, const vector_fp& order) { + m_n[rxn] = k.size(); + switch (k.size()) { + case 1: + m_loc[rxn] = m_c1_list.size(); + m_c1_list.push_back(C1(rxn, k[0], order[0])); + break; + case 2: + m_loc[rxn] = m_c2_list.size(); + m_c2_list.push_back(C2(rxn, k[0], k[1], order[0], order[1])); + break; + case 3: + m_loc[rxn] = m_c3_list.size(); + m_c3_list.push_back(C3(rxn, k[0], k[1], k[2], + order[0], order[1], order[2])); + break; + default: + m_loc[rxn] = m_cn_list.size(); + m_cn_list.push_back(C_AnyN(rxn, k, order)); + } + } - void multiply(const doublereal* input, doublereal* output) const { - _multiply(m_c1_list.begin(), m_c1_list.end(), input, output); - _multiply(m_c2_list.begin(), m_c2_list.end(), input, output); - _multiply(m_c3_list.begin(), m_c3_list.end(), input, output); - _multiply(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void multiply(const doublereal* input, doublereal* output) const { + _multiply(m_c1_list.begin(), m_c1_list.end(), input, output); + _multiply(m_c2_list.begin(), m_c2_list.end(), input, output); + _multiply(m_c3_list.begin(), m_c3_list.end(), input, output); + _multiply(m_cn_list.begin(), m_cn_list.end(), input, output); + } - void power(const doublereal* input, doublereal* output) const { - _power(m_c1_list.begin(), m_c1_list.end(), input, output); - _power(m_c2_list.begin(), m_c2_list.end(), input, output); - _power(m_c3_list.begin(), m_c3_list.end(), input, output); - _power(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void power(const doublereal* input, doublereal* output) const { + _power(m_c1_list.begin(), m_c1_list.end(), input, output); + _power(m_c2_list.begin(), m_c2_list.end(), input, output); + _power(m_c3_list.begin(), m_c3_list.end(), input, output); + _power(m_cn_list.begin(), m_cn_list.end(), input, output); + } - void incrementSpecies(const doublereal* input, doublereal* output) const { - _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); - _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); - _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); - _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void incrementSpecies(const doublereal* input, doublereal* output) const { + _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); + _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); + _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); + _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); + } - void decrementSpecies(const doublereal* input, doublereal* output) const { - _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); - _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); - _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); - _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void decrementSpecies(const doublereal* input, doublereal* output) const { + _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); + _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); + _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); + _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); + } - void incrementReactions(const doublereal* input, doublereal* output) const { - _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); - _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); - _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); - _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void incrementReactions(const doublereal* input, doublereal* output) const { + _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); + _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); + _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); + _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); + } - void decrementReactions(const doublereal* input, doublereal* output) const { - _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); - _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); - _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); - _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); - } + void decrementReactions(const doublereal* input, doublereal* output) const { + _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); + _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); + _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); + _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); + } - private: + private: - vector m_c1_list; - vector m_c2_list; - vector m_c3_list; - vector m_cn_list; - /** - * Mapping with the Reaction Number as key and the Number of species - * as the value. - */ - map m_n; - /** - * Mapping with the Reaction Number as key and the placement in the - * vector of reactions list( i.e., m_c1_list[]) as key - */ - map m_loc; - }; + vector m_c1_list; + vector m_c2_list; + vector m_c3_list; + vector m_cn_list; + /** + * Mapping with the Reaction Number as key and the Number of species + * as the value. + */ + map m_n; + /** + * Mapping with the Reaction Number as key and the placement in the + * vector of reactions list( i.e., m_c1_list[]) as key + */ + map m_loc; + }; - class StoichWriter { - public: + class StoichWriter { + public: - StoichWriter() {} + StoichWriter() {} - void add(int rxn, const vector_int& k) { - int n, nn = k.size(); - for (n = 0; n < nn; n++) { - if (m_mult[rxn] != "") m_mult[rxn] += " * "; - m_mult[rxn] += "c[" + int2str(k[n]) + "]"; - m_is[k[n]] += " + rop[" + int2str(rxn) + "]"; - m_ds[k[n]] += " - rop[" + int2str(rxn) + "]"; - m_ir[rxn] += " + grt[" + int2str(k[n]) + "]"; - m_dr[rxn] += " - grt[" + int2str(k[n]) + "]"; - } - } - - void writeIncSpec(ostream& s, int nsp) { - int k; - for (k = 0; k < nsp; k++) { - s << "out[" << k << "] = " << m_is[k] << ";" << endl; - } - } - - string mult(int rxn) { return m_mult[rxn]; } - string incrSpec(int k) { return m_is[k]; } - string decrSpec(int k) { return m_ds[k]; } - string incrRxn(int rxn) { return m_ir[rxn]; } - string decrRxn(int rxn) { return m_dr[rxn]; } - - private: - map m_mult, m_ir, m_dr, m_is, m_ds; - }; + void add(int rxn, const vector_int& k) { + int n, nn = k.size(); + for (n = 0; n < nn; n++) { + if (m_mult[rxn] != "") m_mult[rxn] += " * "; + m_mult[rxn] += "c[" + int2str(k[n]) + "]"; + m_is[k[n]] += " + rop[" + int2str(rxn) + "]"; + m_ds[k[n]] += " - rop[" + int2str(rxn) + "]"; + m_ir[rxn] += " + grt[" + int2str(k[n]) + "]"; + m_dr[rxn] += " - grt[" + int2str(k[n]) + "]"; + } } + void writeIncSpec(ostream& s, int nsp) { + int k; + for (k = 0; k < nsp; k++) { + s << "out[" << k << "] = " << m_is[k] << ";" << endl; + } + } + + string mult(int rxn) { return m_mult[rxn]; } + string incrSpec(int k) { return m_is[k]; } + string decrSpec(int k) { return m_ds[k]; } + string incrRxn(int rxn) { return m_ir[rxn]; } + string decrRxn(int rxn) { return m_dr[rxn]; } + + private: + map m_mult, m_ir, m_dr, m_is, m_ds; + }; +} + #endif diff --git a/Cantera/src/config.h b/Cantera/src/config.h index b891e9df9..798526551 100755 --- a/Cantera/src/config.h +++ b/Cantera/src/config.h @@ -1,59 +1,3 @@ -/* ../config.h. Generated automatically by configure. */ -// -// Run the 'configure' script to generate 'config.h' from this input file. -// -#ifndef CT_CONFIG_H -#define CT_CONFIG_H +#include "../../config.h" -//------------------------ Fortran settings -------------------// - - -// define types doublereal, integer, and ftnlen to match the -// corresponding Fortran data types on your system. The defaults -// are OK for most systems - -typedef double doublereal; // Fortran double precision -typedef int integer; // Fortran integer -typedef int ftnlen; // Fortran hidden string length type - - -// Fortran compilers pass character strings in argument lists by -// adding a hidden argement with the length of the string. Some -// compilers add the hidden length argument immediately after the -// CHARACTER variable being passed, while others put all of the hidden -// length arguments at the end of the argument list. Define this if -// the lengths are at the end of the argument list. This is usually the -// case for most unix Fortran compilers, but is (by default) false for -// Visual Fortran under Windows. -#define STRING_LEN_AT_END - - -// Define this if Fortran adds a trailing underscore to names in object files. -// For linux and most unix systems, this is the case. -#define FTN_TRAILING_UNDERSCORE - - -//-------- LAPACK / BLAS --------- - -// Define if you are using LAPACK and BLAS from the Intel Math Kernel -// Library -/* #undef HAVE_INTEL_MKL */ - -#define LAPACK_FTN_STRING_LEN_AT_END 1 -#define LAPACK_NAMES_LOWERCASE 1 -#define LAPACK_FTN_TRAILING_UNDERSCORE 1 - - -//--------- Cantera -------------- - - -//--------- CKReader ------------- - - - -//--------- CtLib ---------------- - - - -#endif diff --git a/Cantera/src/ct_defs.h b/Cantera/src/ct_defs.h index 0c9c3e7d8..239322bdd 100755 --- a/Cantera/src/ct_defs.h +++ b/Cantera/src/ct_defs.h @@ -21,7 +21,7 @@ #include "config.h" // STL includes -//#include +#include #include #include #include diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp index 3ad9871d3..daeb9dd49 100755 --- a/Cantera/src/misc.cpp +++ b/Cantera/src/misc.cpp @@ -238,21 +238,7 @@ namespace Cantera { /** - * Set the default directories for input files. Four directories are - * added to the search path used by findInputFile. These are - * 'data', 'data/inputs', 'data/thermo', and - * 'data/transport'. These names are for convenience only - - * findInputFile searches all of them, independent of the type of - * file. The location of the 'data' directory depends on how - * environment variables are set. If CANTERA_DATA_DIR is set, then - * this will be used instead of 'data'. In addition, if - * WIN_CANTERA_ROOT or CANTERA_ROOT are set, then 'data' is - * assumed to be a top-level subdirectory. WIN_CANTERA_ROOT should - * only be set on PCs, and should be in 'DOS' format, for example - * 'C:\CANTERA'. CANTERA_ROOT, on the other hand, should be in - * unix-like format ('/home/usr/cantera'). This allows Cantera to - * be built on PCs using a unix-like environment (Cygwin) and - * compiler (g++), as well as using Win32 compilers. + * Set the default directories for input data files. */ void setDefaultDirectories() { appinit(); @@ -261,15 +247,16 @@ namespace Cantera { // always look in the local directory first dirs.push_back("."); + #ifdef WIN32 - /* - * Under Windows, the Cantera setup utility puts data files in - * a directory 'Cantera\data' below the one the environment - * variable COMMONPROGRAMFILES points to. (This is usually - * C:\Program Files\Common Files.) If this environment - * variable is defined, then this directory is assumed to - * exist and is added to the search path. - */ + // + // Under Windows, the Cantera setup utility puts data files in + // a directory 'Cantera\data' below the one the environment + // variable COMMONPROGRAMFILES points to. (This is usually + // C:\Program Files\Common Files.) If this environment + // variable is defined, then this directory is assumed to + // exist and is added to the search path. + // const char* comfiles = getenv("COMMONPROGRAMFILES"); if (comfiles != 0) { string cfiles = string(comfiles); @@ -285,6 +272,18 @@ namespace Cantera { dirs.push_back(tmpldir); } #endif + +#ifdef DARWIN + // + // add a default data location for Mac OS X + // + dirs.push_back("/Applications/Cantera/Data"); +#endif + + // + // if environment variable CANTERA_DATA is defined, then add + // it to the search path + // if (getenv("CANTERA_DATA") != 0) { string datadir = string(getenv("CANTERA_DATA")); dirs.push_back(datadir); @@ -299,6 +298,7 @@ namespace Cantera { string datadir = string(CANTERA_ROOT) + "/data"; dirs.push_back(datadir); #endif + } diff --git a/Cantera/src/units.h b/Cantera/src/units.h index 1ddd437a1..c32029408 100644 --- a/Cantera/src/units.h +++ b/Cantera/src/units.h @@ -8,6 +8,7 @@ #define CT_UNITS_H #include "ct_defs.h" +#include "ctexceptions.h" namespace Cantera {