From 70fb26459e146a1a2759744a906f2c70d21cb74c Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 16 Dec 2008 20:32:18 +0000 Subject: [PATCH] Fleshed out the type() and ID() member functions. started adding a duplMyselfAsKinetics() capability. This is unfinished. Doxygen updates --- Cantera/src/kinetics/EdgeKinetics.h | 4 +- Cantera/src/kinetics/GRI_30_Kinetics.h | 3 +- Cantera/src/kinetics/GasKinetics.h | 3 +- Cantera/src/kinetics/InterfaceKinetics.cpp | 88 + Cantera/src/kinetics/InterfaceKinetics.h | 37 +- Cantera/src/kinetics/Kinetics.cpp | 434 ++--- Cantera/src/kinetics/Kinetics.h | 1727 ++++++++++---------- 7 files changed, 1261 insertions(+), 1035 deletions(-) diff --git a/Cantera/src/kinetics/EdgeKinetics.h b/Cantera/src/kinetics/EdgeKinetics.h index 00bb6d15f..b19e7b3e6 100644 --- a/Cantera/src/kinetics/EdgeKinetics.h +++ b/Cantera/src/kinetics/EdgeKinetics.h @@ -41,13 +41,13 @@ namespace Cantera { * Identifies the subclass of the Kinetics manager type. * These are listed in mix_defs.h. */ - virtual int ID() { return cEdgeKinetics; } + virtual int ID() const { return cEdgeKinetics; } /** * Identifies the subclass of the Kinetics manager type. * These are listed in mix_defs.h. */ - virtual int type() { return cEdgeKinetics; } + virtual int type() const { return cEdgeKinetics; } // defined in InterfaceKinetics.cpp virtual void finalize(); diff --git a/Cantera/src/kinetics/GRI_30_Kinetics.h b/Cantera/src/kinetics/GRI_30_Kinetics.h index 0d6d16fca..094c61289 100755 --- a/Cantera/src/kinetics/GRI_30_Kinetics.h +++ b/Cantera/src/kinetics/GRI_30_Kinetics.h @@ -30,7 +30,8 @@ namespace Cantera { /// Destructor. virtual ~GRI_30_Kinetics(){} - virtual int ID() { return cGRI_30_Kinetics; } + virtual int ID() const { return cGRI_30_Kinetics; } + virtual int type() const { return cGRI_30_Kinetics; } virtual void getNetProductionRates(doublereal* net) { gri30_updateROP(); diff --git a/Cantera/src/kinetics/GasKinetics.h b/Cantera/src/kinetics/GasKinetics.h index b8c3324ef..2e74ea96d 100755 --- a/Cantera/src/kinetics/GasKinetics.h +++ b/Cantera/src/kinetics/GasKinetics.h @@ -89,7 +89,8 @@ namespace Cantera { /// Destructor. virtual ~GasKinetics(); - virtual int ID() { return cGasKinetics; } + virtual int ID() const { return cGasKinetics; } + virtual int type() const { return cGasKinetics; } virtual doublereal reactantStoichCoeff(int k, int i) const { return m_rrxn[k][i]; diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index 80f1c353f..c5fd96a2f 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -65,6 +65,94 @@ namespace Cantera { } } + // Copy Constructor for the %InterfaceKinetics object. + /* + * Currently, this is not fully implemented. If called it will + * throw an exception. + */ + InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics &right) : + Kinetics(), + m_kk(0), + m_redo_rates(false), + m_nirrev(0), + m_nrev(0), + m_surf(0), + m_integrator(0), + m_finalized(false), + m_has_coverage_dependence(false), + m_has_electrochem_rxns(false), + m_ioFlag(0) + { + m_kdata = new InterfaceKineticsData; + m_kdata->m_temp = 0.0; + /* + * Call the assignment operator + */ + *this = operator=(right); + } + + // Assignment operator + /* + * This is NOT a virtual function. + * + * @param right Reference to %Kinetics object to be copied into the + * current one. + */ + InterfaceKinetics& InterfaceKinetics:: + operator=(const InterfaceKinetics &right) { + /* + * Check for self assignment. + */ + if (this == &right) return *this; + + Kinetics::operator=(right); + + m_kk = right.m_kk; + m_revindex = right.m_revindex; + m_rates = right.m_rates; + m_redo_rates = right.m_redo_rates; + m_index = right.m_index; + m_irrev = right.m_irrev; + m_rxnstoich = right.m_rxnstoich; + m_nirrev = right.m_nirrev; + m_nrev = right.m_nrev; + m_rrxn = right.m_rrxn; + m_prxn = right.m_prxn; + m_rxneqn = right.m_rxneqn; + *m_kdata = *right.m_kdata; // needs to be developed + m_mu0 = right.m_mu0; + m_phi = right.m_phi; + m_pot = right.m_pot; + m_rwork = right.m_rwork; + m_E = right.m_E; + m_surf = right.m_surf; //DANGER - shallow copy + m_integrator = right.m_integrator; //DANGER - shallow copy + m_beta = right.m_beta; + m_ctrxn = right.m_ctrxn; + m_finalized = right.m_finalized; + m_has_coverage_dependence = right.m_has_coverage_dependence; + m_has_electrochem_rxns = right.m_has_electrochem_rxns; + m_ioFlag = right.m_ioFlag; + + return *this; + } + + + // Duplication routine for objects which inherit from + // Kinetics + /* + * This virtual routine can be used to duplicate %Kinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + Kinetics *InterfaceKinetics::duplMyselfAsKinetics() const { + InterfaceKinetics* tp = new InterfaceKinetics(*this); + return dynamic_cast(tp); + } + /** * Update properties that depend on temperature diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 1e3dd76c2..ffc7259b7 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -88,8 +88,38 @@ namespace Cantera { /// Destructor. virtual ~InterfaceKinetics(); - virtual int ID() { return cInterfaceKinetics; } - virtual int type() { return cInterfaceKinetics; } + //! Copy Constructor for the %Kinetics object. + /*! + * Currently, this is not fully implemented. If called it will + * throw an exception. + */ + InterfaceKinetics(const InterfaceKinetics &right); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %Kinetics object to be copied into the + * current one. + */ + InterfaceKinetics& operator=(const InterfaceKinetics &right); + + + //! Duplication routine for objects which inherit from + //! Kinetics + /*! + * This virtual routine can be used to duplicate %InterfaceKinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + virtual Kinetics *duplMyselfAsKinetics() const; + + + virtual int ID() const { return cInterfaceKinetics; } + virtual int type() const { return cInterfaceKinetics; } /** * Set the electric potential in the nth phase @@ -487,6 +517,7 @@ namespace Cantera { * product stoichiometric coefficient for the species being the value. */ mutable std::vector > m_prxn; + //! String expression for each rxn /*! * Vector of strings of length m_ii, the number of @@ -525,7 +556,7 @@ namespace Cantera { * for all of the species in the kinetics object * * Length = m_k - * units = J/kmol + * units = J/kmol */ vector_fp m_mu0; diff --git a/Cantera/src/kinetics/Kinetics.cpp b/Cantera/src/kinetics/Kinetics.cpp index 3ac5b17ec..38ea75fe0 100644 --- a/Cantera/src/kinetics/Kinetics.cpp +++ b/Cantera/src/kinetics/Kinetics.cpp @@ -6,11 +6,15 @@ * * Kinetics managers calculate rates of progress of species due to homogeneous or heterogeneous kinetics. */ +/* + * $Date$ + * $Revision$ + */ // Copyright 2001-2004 California Institute of Technology - + #include "InterfaceKinetics.h" #include "SurfPhase.h" #include "StoichManager.h" @@ -25,197 +29,261 @@ using namespace std; namespace Cantera { - Kinetics::Kinetics() : m_ii(0), m_thermo(0), - m_index(-1), m_surfphase(-1), m_rxnphase(-1), - m_mindim(4) {} + Kinetics::Kinetics() : m_ii(0), m_thermo(0), + m_index(-1), m_surfphase(-1), m_rxnphase(-1), + m_mindim(4) {} + Kinetics::~Kinetics(){} + + + // Copy Constructor for the %Kinetics object. + /* + * Currently, this is not fully implemented. If called it will + * throw an exception. + */ + Kinetics::Kinetics(const Kinetics &right) : + m_ii(0), + m_thermo(0), + m_index(-1), + m_surfphase(-1), + m_rxnphase(-1), + m_mindim(4) + { /* - Kinetics::Kinetics(thermo_t* thermo) : - m_ii(0), - m_index(-1), - m_surfphase(-1) , - m_rxnphase(0), m_mindim(4) - { - if (thermo) { - addPhase(*thermo); - } - deprecatedMethod("Kinetics","Kinetics(thermo_t*)","Kinetics()"); - removeAtVersion("Kinetics(thermo_t*)","1.6.0"); - } - */ - - Kinetics::~Kinetics(){} - - /** - * Takes as input an array of properties for all species in the - * mechanism and copies those values beloning to a particular - * phase to the output array. - * @param data Input data array. - * @param phase Pointer to one of the phase objects participating - * in this reaction mechanism - * @param phase_data Output array where the values for the the - * specified phase are to be written. + * Call the assignment operator */ - void Kinetics::selectPhase(const doublereal* data, const thermo_t* phase, - doublereal* phase_data) { - int n, nsp, np = nPhases(); - for (n = 0; n < np; n++) { - if (phase == m_thermo[n]) { - nsp = phase->nSpecies(); - copy(data + m_start[n], - data + m_start[n] + nsp, phase_data); - return; - } - } - throw CanteraError("Kinetics::selectPhase", "Phase not found."); - } - - - /** - * kineticsSpeciesName(): - * - * Return the string name of the kth species in the kinetics - * manager. k is an integer from 0 to ktot - 1, where ktot is - * the number of species in the kinetics manager, which is the - * sum of the number of species in all phases participating in - * the kinetics manager. If k is out of bounds, the string - * "" is returned. + *this = operator=(right); + } + + // Assignment operator + /* + * This is NOT a virtual function. + * + * @param right Reference to %Kinetics object to be copied into the + * current one. + */ + Kinetics& Kinetics:: + operator=(const Kinetics &right) { + /* + * Check for self assignment. */ - string Kinetics::kineticsSpeciesName(int k) const { - int np = m_start.size(); - for (int n = np-1; n >= 0; n--) { - if (k >= m_start[n]) { - return thermo(n).speciesName(k - m_start[n]); - } - } - return ""; + if (this == &right) return *this; + + m_ii = right.m_ii; + m_perturb = right.m_perturb; + m_reactants = right.m_reactants; + m_products = right.m_products; + + m_thermo = right.m_thermo; // DANGER -> shallow pointer copy + + m_start = right.m_start; + m_phaseindex = right.m_phaseindex; + m_index = right.m_index; + m_surfphase = right.m_surfphase; + m_rxnphase = right.m_rxnphase; + m_mindim = right.m_mindim; + m_dummygroups = right.m_dummygroups; + + return *this; + } + + + // Duplication routine for objects which inherit from + // Kinetics + /* + * This virtual routine can be used to duplicate %Kinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + Kinetics *Kinetics::duplMyselfAsKinetics() const { + Kinetics* tp = new Kinetics(*this); + return tp; + } + + + + int Kinetics::ID() const { + return 0; + } + + int Kinetics::type() const { + return 0; + } + + + /** + * Takes as input an array of properties for all species in the + * mechanism and copies those values beloning to a particular + * phase to the output array. + * @param data Input data array. + * @param phase Pointer to one of the phase objects participating + * in this reaction mechanism + * @param phase_data Output array where the values for the the + * specified phase are to be written. + */ + void Kinetics::selectPhase(const doublereal* data, const thermo_t* phase, + doublereal* phase_data) { + int n, nsp, np = nPhases(); + for (n = 0; n < np; n++) { + if (phase == m_thermo[n]) { + nsp = phase->nSpecies(); + copy(data + m_start[n], + data + m_start[n] + nsp, phase_data); + return; + } + } + throw CanteraError("Kinetics::selectPhase", "Phase not found."); + } + + + /** + * kineticsSpeciesName(): + * + * Return the string name of the kth species in the kinetics + * manager. k is an integer from 0 to ktot - 1, where ktot is + * the number of species in the kinetics manager, which is the + * sum of the number of species in all phases participating in + * the kinetics manager. If k is out of bounds, the string + * "" is returned. + */ + string Kinetics::kineticsSpeciesName(int k) const { + int np = m_start.size(); + for (int n = np-1; n >= 0; n--) { + if (k >= m_start[n]) { + return thermo(n).speciesName(k - m_start[n]); + } + } + return ""; + } + + /** + * kineticsSpeciesIndex(): + * + * This routine will look up a species number based on + * the input string nm. The lookup of species will + * occur for all phases listed in the kinetics object, + * unless the string ph refers to a specific phase of + * the object. + * + * return + * - If a match is found, the position in the species list + * is returned. + * - If a specific phase is specified and no match is found, + * the value -1 is returned. + * - If no match is found in any phase, the value -2 is returned. + */ + int Kinetics::kineticsSpeciesIndex(std::string nm, std::string ph) const { + int np = static_cast(m_thermo.size()); + int k; + string id; + for (int n = 0; n < np; n++) { + id = thermo(n).id(); + if (ph == id) { + k = thermo(n).speciesIndex(nm); + if (k < 0) return -1; + return k + m_start[n]; + } + else if (ph == "") { + /* + * Call the speciesIndex() member function of the + * ThermoPhase object to find a match. + */ + k = thermo(n).speciesIndex(nm); + if (k >= 0) return k + m_start[n]; + } + } + return -2; + } + + /** + * This function looks up the string name of a species and + * returns a reference to the ThermoPhase object of the + * phase where the species resides. + * Will throw an error if the species string doesn't match. + */ + thermo_t& Kinetics::speciesPhase(std::string nm) { + int np = static_cast(m_thermo.size()); + int k; + string id; + for (int n = 0; n < np; n++) { + k = thermo(n).speciesIndex(nm); + if (k >= 0) return thermo(n); + } + throw CanteraError("speciesPhase", "unknown species "+nm); + } + + /** + * This function takes as an argument the kineticsSpecies index + * (i.e., the list index in the list of species in the kinetics + * manager) and returns the index of the phase owning the + * species. + */ + int Kinetics::speciesPhaseIndex(int k) { + int np = m_start.size(); + for (int n = np-1; n >= 0; n--) { + if (k >= m_start[n]) { + return n; + } + } + throw CanteraError("speciesPhaseIndex", + "illegal species index: "+int2str(k)); + } + + /* + * Add a phase to the kinetics manager object. This must + * be done before the function init() is called or + * before any reactions are input. + * The following fields are updated: + * m_start -> vector of integers, containing the + * starting position of the species for + * each phase in the kinetics mechanism. + * m_surfphase -> index of the surface phase. + * m_thermo -> vector of pointers to ThermoPhase phases + * that participate in the kinetics + * mechanism. + * m_phaseindex -> map containing the string id of each + * ThermoPhase phase as a key and the + * index of the phase within the kinetics + * manager object as the value. + */ + void Kinetics::addPhase(thermo_t& thermo) { + + // if not the first thermo object, set the start position + // to that of the last object added + the number of its species + if (m_thermo.size() > 0) { + m_start.push_back(m_start.back() + + m_thermo.back()->nSpecies()); + } + // otherwise start at 0 + else { + m_start.push_back(0); } - /** - * kineticsSpeciesIndex(): - * - * This routine will look up a species number based on - * the input string nm. The lookup of species will - * occur for all phases listed in the kinetics object, - * unless the string ph refers to a specific phase of - * the object. - * - * return - * - If a match is found, the position in the species list - * is returned. - * - If a specific phase is specified and no match is found, - * the value -1 is returned. - * - If no match is found in any phase, the value -2 is returned. - */ - int Kinetics::kineticsSpeciesIndex(std::string nm, std::string ph) const { - int np = static_cast(m_thermo.size()); - int k; - string id; - for (int n = 0; n < np; n++) { - id = thermo(n).id(); - if (ph == id) { - k = thermo(n).speciesIndex(nm); - if (k < 0) return -1; - return k + m_start[n]; - } - else if (ph == "") { - /* - * Call the speciesIndex() member function of the - * ThermoPhase object to find a match. - */ - k = thermo(n).speciesIndex(nm); - if (k >= 0) return k + m_start[n]; - } - } - return -2; + // the phase with lowest dimensionality is assumed to be the + // phase/interface at which reactions take place + if (thermo.nDim() <= m_mindim) { + m_mindim = thermo.nDim(); + m_rxnphase = nPhases(); } - /** - * This function looks up the string name of a species and - * returns a reference to the ThermoPhase object of the - * phase where the species resides. - * Will throw an error if the species string doesn't match. - */ - thermo_t& Kinetics::speciesPhase(std::string nm) { - int np = static_cast(m_thermo.size()); - int k; - string id; - for (int n = 0; n < np; n++) { - k = thermo(n).speciesIndex(nm); - if (k >= 0) return thermo(n); - } - throw CanteraError("speciesPhase", "unknown species "+nm); - } - - /** - * This function takes as an argument the kineticsSpecies index - * (i.e., the list index in the list of species in the kinetics - * manager) and returns the index of the phase owning the - * species. - */ - int Kinetics::speciesPhaseIndex(int k) { - int np = m_start.size(); - for (int n = np-1; n >= 0; n--) { - if (k >= m_start[n]) { - return n; - } - } - throw CanteraError("speciesPhaseIndex", - "illegal species index: "+int2str(k)); - } - - /** - * Add a phase to the kinetics manager object. This must - * be done before the function init() is called or - * before any reactions are input. - * The following fields are updated: - * m_start -> vector of integers, containing the - * starting position of the species for - * each phase in the kinetics mechanism. - * m_surfphase -> index of the surface phase. - * m_thermo -> vector of pointers to ThermoPhase phases - * that participate in the kinetics - * mechanism. - * m_phaseindex -> map containing the string id of each - * ThermoPhase phase as a key and the - * index of the phase within the kinetics - * manager object as the value. - */ - void Kinetics::addPhase(thermo_t& thermo) { - - // if not the first thermo object, set the start position - // to that of the last object added + the number of its species - if (m_thermo.size() > 0) { - m_start.push_back(m_start.back() - + m_thermo.back()->nSpecies()); - } - // otherwise start at 0 - else { - m_start.push_back(0); - } - - // the phase with lowest dimensionality is assumed to be the - // phase/interface at which reactions take place - if (thermo.nDim() <= m_mindim) { - m_mindim = thermo.nDim(); - m_rxnphase = nPhases(); - } - - // there should only be one surface phase - int ptype = -100; - if (type() == cEdgeKinetics) ptype = cEdge; - else if (type() == cInterfaceKinetics) ptype = cSurf; - if (thermo.eosType() == ptype) { - // if (m_surfphase >= 0) { - // throw CanteraError("Kinetics::addPhase", - // "cannot add more than one surface phase"); - // } - m_surfphase = nPhases(); - m_rxnphase = nPhases(); - } - m_thermo.push_back(&thermo); - m_phaseindex[m_thermo.back()->id()] = nPhases(); + // there should only be one surface phase + int ptype = -100; + if (type() == cEdgeKinetics) ptype = cEdge; + else if (type() == cInterfaceKinetics) ptype = cSurf; + if (thermo.eosType() == ptype) { + // if (m_surfphase >= 0) { + // throw CanteraError("Kinetics::addPhase", + // "cannot add more than one surface phase"); + // } + m_surfphase = nPhases(); + m_rxnphase = nPhases(); } + m_thermo.push_back(&thermo); + m_phaseindex[m_thermo.back()->id()] = nPhases(); + } //! Private function of the class Kinetics, indicating that a function diff --git a/Cantera/src/kinetics/Kinetics.h b/Cantera/src/kinetics/Kinetics.h index e849984f4..25ca5305f 100755 --- a/Cantera/src/kinetics/Kinetics.h +++ b/Cantera/src/kinetics/Kinetics.h @@ -4,7 +4,6 @@ * module documentation (see \ref kineticsmgr and class * \link Cantera::Kinetics Kinetics\endlink). * - * $Author$ * $Date$ * $Revision$ */ @@ -20,951 +19,989 @@ namespace Cantera { - // forward references - class ReactionData; + // forward references + class ReactionData; + + /** + * @defgroup chemkinetics Chemical Kinetics + */ + + /// @defgroup kineticsmgr Kinetics Managers + /// @section kinmodman Models and Managers + /// + /// A kinetics manager is a C++ class that implements a kinetics + /// model; a kinetics model is a set of mathematical equation + /// describing how various kinetic quanities are to be computed -- + /// reaction rates, species production rates, etc. Many different + /// kinetics models might be defined to handle different types of + /// kinetic processes. For example, one kinetics model might use + /// expressions valid for elementary reactions in ideal gas + /// mixtures. It might, for example, require the reaction orders + /// to be integral and equal to the forward stoichiometric + /// coefficients, require that each reaction be reversible with a + /// reverse rate satisfying detailed balance, include + /// pressure-dependent unimolecular reactions, etc. Another + /// kinetics model might be designed for heterogeneous chemistry + /// at interfaces, and might allow empirical reaction orders, + /// coverage-dependent activation energies, irreversible + /// reactions, and include effects of potential differences across + /// the interface on reaction rates. + /// + /// A kinetics manager implements a kinetics model. Since the + /// model equations may be complex and expensive to evaluate, a + /// kinetics manager may adopt various strategies to 'manage' the + /// computation and evaluate the expressions efficiently. For + /// example, if there are rate coefficients or other quantities + /// that depend only on temperature, a manager class may choose to + /// store these quantities internally, and re-evaluate them only + /// when the temperature has actually changed. Or a manager + /// designed for use with reaction mechanisms with a few repeated + /// activation energies might precompute the terms \f$ exp(-E/RT) + /// \f$, instead of evaluating the exponential repeatedly for each + /// reaction. There are many other possible 'management styles', + /// each of which might be better suited to some reaction + /// mechanisms than others. + /// + /// But however a manager structures the internal computation, the + /// tasks the manager class must perform are, for the most part, + /// the same. It must be able to compute reaction rates, species + /// production rates, equilibrium constants, etc. Therefore, all + /// kinetics manager classes should have a common set of public + /// methods, but differ in how they implement these methods. + /// + /// A kinetics manager computes reaction rates of progress, + /// species production rates, equilibrium constants, and similar + /// quantities for a reaction mechanism. All kinetics manager + /// classes derive from class Kinetics, which defines a common + /// public interface for all kinetics managers. Each derived class + /// overloads the virtual methods of Kinetics to implement a + /// particular kinetics model. + /// + /// For example, class GasKinetics implements reaction rate + /// expressions appropriate for homogeneous reactions in ideal gas + /// mixtures, and class InterfaceKinetics implements expressions + /// appropriate for heterogeneous mechanisms at interfaces, + /// including how to handle reactions involving charged species of + /// phases with different electric potentials --- something that + /// class GasKinetics doesn't deal with at all. + /// + /// Kinetics managers may be also created that hard-wire a + /// particular reaction mechanism in C++ code. This can often + /// result in faster performance. An example of this is the + /// kinetics manager GRI30_Kinetics that hard-wires the rate + /// expressions for the natural gas combustion mechanism GRI-3.0. + /// + /// Many of the methods of class Kinetics write into arrays the + /// values of some quantity for each species, for example the net + /// production rate. These methods always write the results into + /// flat arrays, ordered by phase in the order the phase was + /// added, and within a phase in the order the species were added + /// to the phase (which is the same ordering as in the input + /// file). Example: suppose a heterogeneous mechanism involves + /// three phases -- a bulk phase 'a', another bulk phase 'b', and + /// the surface phase 'a:b' at the a/b interface. Phase 'a' + /// contains 12 species, phase 'b' contains 3, and at the + /// interface there are 5 adsorbed species defined in phase + /// 'a:b'. Then methods like getNetProductionRates(doublereal* net) + /// will write and output array of length 20, beginning at the location + /// pointed to by 'net'. The first 12 values will be the net production + /// rates for all 12 species of phase 'a' (even if some do not participate + /// in the reactions), the next 3 will be for phase 'b', and finally the + /// net production rates for the surface species will occupy the last + /// 5 locations. + /// @ingroup chemkinetics + + + //! Public interface for kinetics managers. + /*! + * This class serves as a base class to derive 'kinetics + * managers', which are classes that manage homogeneous chemistry + * within one phase, or heterogeneous chemistry at one + * interface. The virtual methods of this class are meant to be + * overloaded in subclasses. The non-virtual methods perform + * generic functions and are implemented in Kinetics. They should + * not be overloaded. Only those methods required by a subclass + * need to be overloaded; the rest will throw exceptions if + * called. + * + * When the nomenclature "kinetics species index" is used below, + * this means that the species index ranges over all species in + * all phases handled by the kinetics manager. + * + * @ingroup kineticsmgr + */ + class Kinetics { + + public: + + //! typedef for ThermoPhase + typedef ThermoPhase thermo_t; /** - * @defgroup chemkinetics Chemical Kinetics + * @name Constructors and General Information about Mechanism */ + //@{ - /// @defgroup kineticsmgr Kinetics Managers - /// @section kinmodman Models and Managers - /// - /// A kinetics manager is a C++ class that implements a kinetics - /// model; a kinetics model is a set of mathematical equation - /// describing how various kinetic quanities are to be computed -- - /// reaction rates, species production rates, etc. Many different - /// kinetics models might be defined to handle different types of - /// kinetic processes. For example, one kinetics model might use - /// expressions valid for elementary reactions in ideal gas - /// mixtures. It might, for example, require the reaction orders - /// to be integral and equal to the forward stoichiometric - /// coefficients, require that each reaction be reversible with a - /// reverse rate satisfying detailed balance, include - /// pressure-dependent unimolecular reactions, etc. Another - /// kinetics model might be designed for heterogeneous chemistry - /// at interfaces, and might allow empirical reaction orders, - /// coverage-dependent activation energies, irreversible - /// reactions, and include effects of potential differences across - /// the interface on reaction rates. - /// - /// A kinetics manager implements a kinetics model. Since the - /// model equations may be complex and expensive to evaluate, a - /// kinetics manager may adopt various strategies to 'manage' the - /// computation and evaluate the expressions efficiently. For - /// example, if there are rate coefficients or other quantities - /// that depend only on temperature, a manager class may choose to - /// store these quantities internally, and re-evaluate them only - /// when the temperature has actually changed. Or a manager - /// designed for use with reaction mechanisms with a few repeated - /// activation energies might precompute the terms \f$ exp(-E/RT) - /// \f$, instead of evaluating the exponential repeatedly for each - /// reaction. There are many other possible 'management styles', - /// each of which might be better suited to some reaction - /// mechanisms than others. - /// - /// But however a manager structures the internal computation, the - /// tasks the manager class must perform are, for the most part, - /// the same. It must be able to compute reaction rates, species - /// production rates, equilibrium constants, etc. Therefore, all - /// kinetics manager classes should have a common set of public - /// methods, but differ in how they implement these methods. - /// - /// A kinetics manager computes reaction rates of progress, - /// species production rates, equilibrium constants, and similar - /// quantities for a reaction mechanism. All kinetics manager - /// classes derive from class Kinetics, which defines a common - /// public interface for all kinetics managers. Each derived class - /// overloads the virtual methods of Kinetics to implement a - /// particular kinetics model. - /// - /// For example, class GasKinetics implements reaction rate - /// expressions appropriate for homogeneous reactions in ideal gas - /// mixtures, and class InterfaceKinetics implements expressions - /// appropriate for heterogeneous mechanisms at interfaces, - /// including how to handle reactions involving charged species of - /// phases with different electric potentials --- something that - /// class GasKinetics doesn't deal with at all. - /// - /// Kinetics managers may be also created that hard-wire a - /// particular reaction mechanism in C++ code. This can often - /// result in faster performance. An example of this is the - /// kinetics manager GRI30_Kinetics that hard-wires the rate - /// expressions for the natural gas combustion mechanism GRI-3.0. - /// - /// Many of the methods of class Kinetics write into arrays the - /// values of some quantity for each species, for example the net - /// production rate. These methods always write the results into - /// flat arrays, ordered by phase in the order the phase was - /// added, and within a phase in the order the species were added - /// to the phase (which is the same ordering as in the input - /// file). Example: suppose a heterogeneous mechanism involves - /// three phases -- a bulk phase 'a', another bulk phase 'b', and - /// the surface phase 'a:b' at the a/b interface. Phase 'a' - /// contains 12 species, phase 'b' contains 3, and at the - /// interface there are 5 adsorbed species defined in phase - /// 'a:b'. Then methods like getNetProductionRates(doublereal* net) - /// will write and output array of length 20, beginning at the location - /// pointed to by 'net'. The first 12 values will be the net production - /// rates for all 12 species of phase 'a' (even if some do not participate - /// in the reactions), the next 3 will be for phase 'b', and finally the - /// net production rates for the surface species will occupy the last - /// 5 locations. - /// @ingroup chemkinetics + /// Default constructor. + Kinetics(); + /// This constructor initializes with a starting phase. + /// @deprecated + // Kinetics(thermo_t* thermo); - //! Public interface for kinetics managers. - /*! - * This class serves as a base class to derive 'kinetics - * managers', which are classes that manage homogeneous chemistry - * within one phase, or heterogeneous chemistry at one - * interface. The virtual methods of this class are meant to be - * overloaded in subclasses. The non-virtual methods perform - * generic functions and are implemented in Kinetics. They should - * not be overloaded. Only those methods required by a subclass - * need to be overloaded; the rest will throw exceptions if - * called. - * - * When the nomenclature "kinetics species index" is used below, - * this means that the species index ranges over all species in - * all phases handled by the kinetics manager. - * - * @ingroup kineticsmgr + /// Destructor. + virtual ~Kinetics(); + + //!Copy Constructor for the %Kinetics object. + /*! + * Currently, this is not fully implemented. If called it will + * throw an exception. */ - class Kinetics { - - public: - - //! typedef for ThermoPhase - typedef ThermoPhase thermo_t; - - /** - * @name Constructors and General Information about Mechanism - */ - //@{ - - /// Default constructor. - Kinetics(); - - /// This constructor initializes with a starting phase. - /// @deprecated - // Kinetics(thermo_t* thermo); - - /// Destructor. - virtual ~Kinetics(); - - /// Identifies the kinetics manager type. Each class derived - /// from Kinetics should overload this method to return a - /// unique integer. Standard values are defined in file - /// mix_defs.h. - virtual int type() { return 0; } - - //! Number of reactions in the reaction mechanism. - int nReactions() const {return m_ii;} - - //@} + Kinetics(const Kinetics &); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %Kinetics object to be copied into the + * current one. + */ + Kinetics& operator=(const Kinetics &right); - /** - * @name Information/Lookup Functions about Phases and Species - */ - //@{ + //! Duplication routine for objects which inherit from + //! Kinetics + /*! + * This virtual routine can be used to duplicate %Kinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + virtual Kinetics *duplMyselfAsKinetics() const; - /** - * The number of phases participating in the reaction - * mechanism. For a homogeneous reaction mechanism, this will - * always return 1, but for a heterogeneous mechanism it will - * return the total number of phases in the mechanism. - */ - int nPhases() const { return static_cast(m_thermo.size()); } + + //! Identifies the subclass of the Kinetics manager type. + /*! + * These are listed in mix_defs.h. + */ + virtual int ID() const; - /** - * Return the phase index of a phase in the list of phases - * defined within the object. - * - * @param ph std::string name of the phase - * - * If a -1 is returned, then the phase is not defined in - * the Kinetics object. - */ - int phaseIndex(std::string ph) { - if (m_phaseindex.find(ph) == m_phaseindex.end()) { - return -1; - } - else { - return m_phaseindex[ph] - 1; - } - } + //! Identifies the kinetics manager type. + /*! + * Each class derived from Kinetics should overload this method to + * return a unique integer. Standard values are defined in file + * mix_defs.h. + */ + virtual int type() const; - /** - * This returns the integer index of the phase which has - * ThermoPhase type cSurf. For heterogeneous mechanisms, this - * identifies the one surface phase. For homogeneous - * mechanisms, this reurns -1. - */ - int surfacePhaseIndex() { return m_surfphase; } + //! Number of reactions in the reaction mechanism. + int nReactions() const {return m_ii;} - /** - * Phase where the reactions occur. For heterogeneous - * mechanisms, one of the phases in the list of phases - * represents the 2D interface or 1D edge at which the - * reactions take place. This method returns the index of the - * phase with the smallest spatial dimension (1, 2, or 3) - * among the list of phases. If there is more than one, the - * index of the first one is returned. For homogeneous - * mechanisms, the value 0 is returned. - */ - int reactionPhaseIndex() { return m_rxnphase; } + //@} - /** - * This method returns a reference to the nth ThermoPhase - * object defined in this kinetics mechanism. It is typically - * used so that member functions of the ThermoPhase object may - * be called. For homogeneous mechanisms, there is only one - * object, and this method can be called without an argument - * to access it. - * - * @param n Index of the ThermoPhase being sought. - */ - thermo_t& thermo(int n=0) { return *m_thermo[n]; } - const thermo_t& thermo(int n=0) const { return *m_thermo[n]; } + /** + * @name Information/Lookup Functions about Phases and Species + */ + //@{ - /** - * This method returns a reference to the nth ThermoPhase - * defined in this kinetics mechanism. - * It is typically used so that member functions of the - * ThermoPhase may be called. @deprecated This method is redundant. - * - * @param n Index of the ThermoPhase being sought. - */ - thermo_t& phase(int n=0) { - deprecatedMethod("Kinetics","phase","thermo"); - return *m_thermo[n]; - } - /** - * This method returns a reference to the nth ThermoPhase - * defined in this kinetics mechanism. - * It is typically used so that member functions of the - * ThermoPhase may be called. @deprecated This method is redundant. - * - * @param n Index of the ThermoPhase being sought. - */ - const thermo_t& phase(int n=0) const { - deprecatedMethod("Kinetics","phase","thermo"); - return *m_thermo[n]; - } + /** + * The number of phases participating in the reaction + * mechanism. For a homogeneous reaction mechanism, this will + * always return 1, but for a heterogeneous mechanism it will + * return the total number of phases in the mechanism. + */ + int nPhases() const { return static_cast(m_thermo.size()); } - /** - * The total number of species in all phases participating in - * the kinetics mechanism. This is useful to dimension arrays - * for use in calls to methods that return the species - * production rates, for example. - */ - int nTotalSpecies() const { - int n=0, np; - np = nPhases(); - for (int p = 0; p < np; p++) n += thermo(p).nSpecies(); - return n; - } + /** + * Return the phase index of a phase in the list of phases + * defined within the object. + * + * @param ph std::string name of the phase + * + * If a -1 is returned, then the phase is not defined in + * the Kinetics object. + */ + int phaseIndex(std::string ph) { + if (m_phaseindex.find(ph) == m_phaseindex.end()) { + return -1; + } + else { + return m_phaseindex[ph] - 1; + } + } - /** - * Returns the starting index of the species in the nth phase - * associated with the reaction mechanism. - * - * @param n Return the index of first species in the nth phase - * associated with the reaction mechanism. - */ - int start(int n) { - deprecatedMethod("Kinetics","start","kineticsSpeciesIndex(0,n)"); - return m_start[n]; - } + /** + * This returns the integer index of the phase which has + * ThermoPhase type cSurf. For heterogeneous mechanisms, this + * identifies the one surface phase. For homogeneous + * mechanisms, this reurns -1. + */ + int surfacePhaseIndex() { return m_surfphase; } + + /** + * Phase where the reactions occur. For heterogeneous + * mechanisms, one of the phases in the list of phases + * represents the 2D interface or 1D edge at which the + * reactions take place. This method returns the index of the + * phase with the smallest spatial dimension (1, 2, or 3) + * among the list of phases. If there is more than one, the + * index of the first one is returned. For homogeneous + * mechanisms, the value 0 is returned. + */ + int reactionPhaseIndex() { return m_rxnphase; } - /** - * The location of species k of phase n in species arrays. - * Kinetics manager classes return species production rates in - * flat arrays, with the species of each phases following one - * another, in the order the phases were added. This method - * is useful to find the value for a particular species of a - * particular phase in arrrays returned from methods like - * getCreationRates that return an array of species-specific - * quantities. - * - * Example: suppose a heterogeneous mechanism involves three - * phases. The first contains 12 species, the second 26, and - * the third 3. Then species arrays must have size at least - * 41, and positions 0 - 11 are the values for the species in - * the first phase, positions 12 - 37 are the values for the - * species in the second phase, etc. Then - * kineticsSpeciesIndex(7, 0) = 7, kineticsSpeciesIndex(4, 1) - * = 16, and kineticsSpeciesIndex(2, 2) = 40. - * - * @param k species index - * @param n phase index for the species - */ - int kineticsSpeciesIndex(int k, int n) const { - return m_start[n] + k; - } + /** + * This method returns a reference to the nth ThermoPhase + * object defined in this kinetics mechanism. It is typically + * used so that member functions of the ThermoPhase object may + * be called. For homogeneous mechanisms, there is only one + * object, and this method can be called without an argument + * to access it. + * + * @param n Index of the ThermoPhase being sought. + */ + thermo_t& thermo(int n=0) { return *m_thermo[n]; } + const thermo_t& thermo(int n=0) const { return *m_thermo[n]; } + + /** + * This method returns a reference to the nth ThermoPhase + * defined in this kinetics mechanism. + * It is typically used so that member functions of the + * ThermoPhase may be called. @deprecated This method is redundant. + * + * @param n Index of the ThermoPhase being sought. + */ + thermo_t& phase(int n=0) { + deprecatedMethod("Kinetics","phase","thermo"); + return *m_thermo[n]; + } + /** + * This method returns a reference to the nth ThermoPhase + * defined in this kinetics mechanism. + * It is typically used so that member functions of the + * ThermoPhase may be called. @deprecated This method is redundant. + * + * @param n Index of the ThermoPhase being sought. + */ + const thermo_t& phase(int n=0) const { + deprecatedMethod("Kinetics","phase","thermo"); + return *m_thermo[n]; + } + + /** + * The total number of species in all phases participating in + * the kinetics mechanism. This is useful to dimension arrays + * for use in calls to methods that return the species + * production rates, for example. + */ + int nTotalSpecies() const { + int n=0, np; + np = nPhases(); + for (int p = 0; p < np; p++) n += thermo(p).nSpecies(); + return n; + } + + /** + * Returns the starting index of the species in the nth phase + * associated with the reaction mechanism. + * + * @param n Return the index of first species in the nth phase + * associated with the reaction mechanism. + */ + int start(int n) { + deprecatedMethod("Kinetics","start","kineticsSpeciesIndex(0,n)"); + return m_start[n]; + } + + + /** + * The location of species k of phase n in species arrays. + * Kinetics manager classes return species production rates in + * flat arrays, with the species of each phases following one + * another, in the order the phases were added. This method + * is useful to find the value for a particular species of a + * particular phase in arrrays returned from methods like + * getCreationRates that return an array of species-specific + * quantities. + * + * Example: suppose a heterogeneous mechanism involves three + * phases. The first contains 12 species, the second 26, and + * the third 3. Then species arrays must have size at least + * 41, and positions 0 - 11 are the values for the species in + * the first phase, positions 12 - 37 are the values for the + * species in the second phase, etc. Then + * kineticsSpeciesIndex(7, 0) = 7, kineticsSpeciesIndex(4, 1) + * = 16, and kineticsSpeciesIndex(2, 2) = 40. + * + * @param k species index + * @param n phase index for the species + */ + int kineticsSpeciesIndex(int k, int n) const { + return m_start[n] + k; + } - //! Return the std::string name of the kth species in the kinetics - //! manager. - /*! - * k is an integer from 0 to ktot - 1, where ktot is - * the number of species in the kinetics manager, which is the - * sum of the number of species in all phases participating in - * the kinetics manager. If k is out of bounds, the std::string - * "" is returned. - * - * @param k species index - */ - std::string kineticsSpeciesName(int k) const; + //! Return the std::string name of the kth species in the kinetics + //! manager. + /*! + * k is an integer from 0 to ktot - 1, where ktot is + * the number of species in the kinetics manager, which is the + * sum of the number of species in all phases participating in + * the kinetics manager. If k is out of bounds, the std::string + * "" is returned. + * + * @param k species index + */ + std::string kineticsSpeciesName(int k) const; - /** - * This routine will look up a species number based on - * the input std::string nm. The lookup of species will - * occur for all phases listed in the kinetics object, - * unless the std::string ph refers to a specific phase of - * the object. - * - * return - * - If a match is found, the position in the species list - * is returned. - * - If a specific phase is specified and no match is found, - * the value -1 is returned. - * - If no match is found in any phase, the value -2 is returned. - * - * @param nm Input string name of the species - * @param ph Input string name of the phase. Defaults to "" - */ - int kineticsSpeciesIndex(std::string nm, std::string ph = "") const; + /** + * This routine will look up a species number based on + * the input std::string nm. The lookup of species will + * occur for all phases listed in the kinetics object, + * unless the std::string ph refers to a specific phase of + * the object. + * + * return + * - If a match is found, the position in the species list + * is returned. + * - If a specific phase is specified and no match is found, + * the value -1 is returned. + * - If no match is found in any phase, the value -2 is returned. + * + * @param nm Input string name of the species + * @param ph Input string name of the phase. Defaults to "" + */ + int kineticsSpeciesIndex(std::string nm, std::string ph = "") const; - /** - * This function looks up the std::string name of a species and - * returns a reference to the ThermoPhase object of the - * phase where the species resides. - * Will throw an error if the species std::string doesn't match. - * - * @param nm String containing the name of the species. - */ - thermo_t& speciesPhase(std::string nm); + /** + * This function looks up the std::string name of a species and + * returns a reference to the ThermoPhase object of the + * phase where the species resides. + * Will throw an error if the species std::string doesn't match. + * + * @param nm String containing the name of the species. + */ + thermo_t& speciesPhase(std::string nm); - /** - * This function takes as an argument the kineticsSpecies index - * (i.e., the list index in the list of species in the kinetics - * manager) and returns the species' owning ThermoPhase object. - * - * @param k Species index - */ - thermo_t& speciesPhase(int k) { - return thermo(speciesPhaseIndex(k)); - } + /** + * This function takes as an argument the kineticsSpecies index + * (i.e., the list index in the list of species in the kinetics + * manager) and returns the species' owning ThermoPhase object. + * + * @param k Species index + */ + thermo_t& speciesPhase(int k) { + return thermo(speciesPhaseIndex(k)); + } - /** - * This function takes as an argument the kineticsSpecies index - * (i.e., the list index in the list of species in the kinetics - * manager) and returns the index of the phase owning the - * species. - * - * @param k Species index - */ - int speciesPhaseIndex(int k); + /** + * This function takes as an argument the kineticsSpecies index + * (i.e., the list index in the list of species in the kinetics + * manager) and returns the index of the phase owning the + * species. + * + * @param k Species index + */ + int speciesPhaseIndex(int k); - //@} + //@} - /** - * @name Reaction Rates Of Progress - */ - //@{ + /** + * @name Reaction Rates Of Progress + */ + //@{ - //! Return the forward rates of progress of the reactions - /*! - * Forward rates of progress. Return the forward rates of - * progress in array fwdROP, which must be dimensioned at - * least as large as the total number of reactions. - * - * @param fwdROP Output vector containing forward rates - * of progress of the reactions. Length: m_ii. - */ - virtual void getFwdRatesOfProgress(doublereal* fwdROP) { - err("getFwdRatesOfProgress"); - } + //! Return the forward rates of progress of the reactions + /*! + * Forward rates of progress. Return the forward rates of + * progress in array fwdROP, which must be dimensioned at + * least as large as the total number of reactions. + * + * @param fwdROP Output vector containing forward rates + * of progress of the reactions. Length: m_ii. + */ + virtual void getFwdRatesOfProgress(doublereal* fwdROP) { + err("getFwdRatesOfProgress"); + } - //! Return the Reverse rates of progress of the reactions - /*! - * Return the reverse rates of - * progress in array revROP, which must be dimensioned at - * least as large as the total number of reactions. - * - * @param revROP Output vector containing reverse rates - * of progress of the reactions. Length: m_ii. - */ - virtual void getRevRatesOfProgress(doublereal* revROP) { - err("getRevRatesOfProgress"); - } + //! Return the Reverse rates of progress of the reactions + /*! + * Return the reverse rates of + * progress in array revROP, which must be dimensioned at + * least as large as the total number of reactions. + * + * @param revROP Output vector containing reverse rates + * of progress of the reactions. Length: m_ii. + */ + virtual void getRevRatesOfProgress(doublereal* revROP) { + err("getRevRatesOfProgress"); + } - /** - * Net rates of progress. Return the net (forward - reverse) - * rates of progress in array netROP, which must be - * dimensioned at least as large as the total number of - * reactions. - * - * @param netROP Output vector of the net ROP. Length: m_ii. - */ - virtual void getNetRatesOfProgress(doublereal* netROP) { - err("getNetRatesOfProgress"); - } + /** + * Net rates of progress. Return the net (forward - reverse) + * rates of progress in array netROP, which must be + * dimensioned at least as large as the total number of + * reactions. + * + * @param netROP Output vector of the net ROP. Length: m_ii. + */ + virtual void getNetRatesOfProgress(doublereal* netROP) { + err("getNetRatesOfProgress"); + } - //! Return a vector of Equilibrium constants. - /*! - * Return the equilibrium constants of - * the reactions in concentration units in array kc, which - * must be dimensioned at least as large as the total number - * of reactions. - * - * @param kc Output vector containing the equilibrium constants. - * Length: m_ii. - */ - virtual void getEquilibriumConstants(doublereal* kc) { - err("getEquilibriumConstants"); - } + //! Return a vector of Equilibrium constants. + /*! + * Return the equilibrium constants of + * the reactions in concentration units in array kc, which + * must be dimensioned at least as large as the total number + * of reactions. + * + * @param kc Output vector containing the equilibrium constants. + * Length: m_ii. + */ + virtual void getEquilibriumConstants(doublereal* kc) { + err("getEquilibriumConstants"); + } - /** - * Change in species properties. Given an array of molar species - * property values \f$ z_k, k = 1, \dots, K \f$, return the - * array of reaction values - * \f[ - * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I. - * \f] - * For example, if this method is called with the array of - * standard-state molar Gibbs free energies for the species, - * then the values returned in array \c deltaProperty would be - * the standard-state Gibbs free energies of reaction for each - * reaction. - * - * @param property Input vector of property value. Length: m_kk. - * @param deltaProperty Output vector of deltaRxn. Length: m_ii. - */ - virtual void getReactionDelta(const doublereal* property, - doublereal* deltaProperty) { - err("getReactionDelta"); - } + /** + * Change in species properties. Given an array of molar species + * property values \f$ z_k, k = 1, \dots, K \f$, return the + * array of reaction values + * \f[ + * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I. + * \f] + * For example, if this method is called with the array of + * standard-state molar Gibbs free energies for the species, + * then the values returned in array \c deltaProperty would be + * the standard-state Gibbs free energies of reaction for each + * reaction. + * + * @param property Input vector of property value. Length: m_kk. + * @param deltaProperty Output vector of deltaRxn. Length: m_ii. + */ + virtual void getReactionDelta(const doublereal* property, + doublereal* deltaProperty) { + err("getReactionDelta"); + } - /** - * Return the vector of values for the reaction gibbs free - * energy change. These values depend upon the concentration - * of the solution. - * - * units = J kmol-1 - * - * @param deltaG Output vector of deltaG's for reactions - * Length: m_ii. - */ - virtual void getDeltaGibbs( doublereal* deltaG) { - err("getDeltaGibbs"); - } + /** + * Return the vector of values for the reaction gibbs free + * energy change. These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + * + * @param deltaG Output vector of deltaG's for reactions + * Length: m_ii. + */ + virtual void getDeltaGibbs( doublereal* deltaG) { + err("getDeltaGibbs"); + } - /** - * Return the vector of values for the reactions change in - * enthalpy. These values depend upon the concentration of - * the solution. - * - * units = J kmol-1 - * - * @param deltaH Output vector of deltaH's for reactions - * Length: m_ii. - */ - virtual void getDeltaEnthalpy( doublereal* deltaH) { - err("getDeltaEnthalpy"); - } + /** + * Return the vector of values for the reactions change in + * enthalpy. These values depend upon the concentration of + * the solution. + * + * units = J kmol-1 + * + * @param deltaH Output vector of deltaH's for reactions + * Length: m_ii. + */ + virtual void getDeltaEnthalpy( doublereal* deltaH) { + err("getDeltaEnthalpy"); + } - /** - * Return the vector of values for the reactions change in - * entropy. These values depend upon the concentration of the - * solution. - * - * units = J kmol-1 Kelvin-1 - * - * @param deltaS Output vector of deltaS's for reactions - * Length: m_ii. - */ - virtual void getDeltaEntropy( doublereal* deltaS) { - err("getDeltaEntropy"); - } + /** + * Return the vector of values for the reactions change in + * entropy. These values depend upon the concentration of the + * solution. + * + * units = J kmol-1 Kelvin-1 + * + * @param deltaS Output vector of deltaS's for reactions + * Length: m_ii. + */ + virtual void getDeltaEntropy( doublereal* deltaS) { + err("getDeltaEntropy"); + } - /** - * Return the vector of values for the reaction standard state - * gibbs free energy change. These values don't depend upon - * the concentration of the solution. - * - * units = J kmol-1 - * - * @param deltaG Output vector of ss deltaG's for reactions - * Length: m_ii. - */ - virtual void getDeltaSSGibbs( doublereal* deltaG) { - err("getDeltaSSGibbs"); - } + /** + * Return the vector of values for the reaction standard state + * gibbs free energy change. These values don't depend upon + * the concentration of the solution. + * + * units = J kmol-1 + * + * @param deltaG Output vector of ss deltaG's for reactions + * Length: m_ii. + */ + virtual void getDeltaSSGibbs( doublereal* deltaG) { + err("getDeltaSSGibbs"); + } - /** - * Return the vector of values for the change in the standard - * state enthalpies of reaction. These values don't depend - * upon the concentration of the solution. - * - * units = J kmol-1 - * - * @param deltaH Output vector of ss deltaH's for reactions - * Length: m_ii. - */ - virtual void getDeltaSSEnthalpy( doublereal* deltaH) { - err("getDeltaSSEnthalpy"); - } + /** + * Return the vector of values for the change in the standard + * state enthalpies of reaction. These values don't depend + * upon the concentration of the solution. + * + * units = J kmol-1 + * + * @param deltaH Output vector of ss deltaH's for reactions + * Length: m_ii. + */ + virtual void getDeltaSSEnthalpy( doublereal* deltaH) { + err("getDeltaSSEnthalpy"); + } - /** - * Return the vector of values for the change in the standard - * state entropies for each reaction. These values don't - * depend upon the concentration of the solution. - * - * units = J kmol-1 Kelvin-1 - * - * @param deltaS Output vector of ss deltaS's for reactions - * Length: m_ii. - */ - virtual void getDeltaSSEntropy( doublereal* deltaS) { - err("getDeltaSSEntropy"); - } + /** + * Return the vector of values for the change in the standard + * state entropies for each reaction. These values don't + * depend upon the concentration of the solution. + * + * units = J kmol-1 Kelvin-1 + * + * @param deltaS Output vector of ss deltaS's for reactions + * Length: m_ii. + */ + virtual void getDeltaSSEntropy( doublereal* deltaS) { + err("getDeltaSSEntropy"); + } - //@} - /** - * @name Species Production Rates - */ - //@{ + //@} + /** + * @name Species Production Rates + */ + //@{ - /** - * Species creation rates [kmol/m^3/s or kmol/m^2/s]. Return the - * species creation rates in array cdot, which must be - * dimensioned at least as large as the total number of - * species in all phases. @see nTotalSpecies. - * - * @param cdot Output vector of creation rates. - * Length: m_kk. - */ - virtual void getCreationRates(doublereal* cdot) { - err("getCreationRates"); - } + /** + * Species creation rates [kmol/m^3/s or kmol/m^2/s]. Return the + * species creation rates in array cdot, which must be + * dimensioned at least as large as the total number of + * species in all phases. @see nTotalSpecies. + * + * @param cdot Output vector of creation rates. + * Length: m_kk. + */ + virtual void getCreationRates(doublereal* cdot) { + err("getCreationRates"); + } - /** - * Species destruction rates [kmol/m^3/s or kmol/m^2/s]. Return - * the species destruction rates in array ddot, which must be - * dimensioned at least as large as the total number of - * species. @see nTotalSpecies. - * - * @param ddot Output vector of destruction rates. - * Length: m_kk. - */ - virtual void getDestructionRates(doublereal* ddot) { - err("getDestructionRates"); - } + /** + * Species destruction rates [kmol/m^3/s or kmol/m^2/s]. Return + * the species destruction rates in array ddot, which must be + * dimensioned at least as large as the total number of + * species. @see nTotalSpecies. + * + * @param ddot Output vector of destruction rates. + * Length: m_kk. + */ + virtual void getDestructionRates(doublereal* ddot) { + err("getDestructionRates"); + } - /** - * Species net production rates [kmol/m^3/s or kmol/m^2/s]. Return - * the species net production rates (creation - destruction) - * in array wdot, which must be dimensioned at least as large - * as the total number of species. @see nTotalSpecies. - * - * @param wdot Output vector of net production rates. - * Length: m_kk. - */ - virtual void getNetProductionRates(doublereal* wdot) { - err("getNetProductionRates"); - } + /** + * Species net production rates [kmol/m^3/s or kmol/m^2/s]. Return + * the species net production rates (creation - destruction) + * in array wdot, which must be dimensioned at least as large + * as the total number of species. @see nTotalSpecies. + * + * @param wdot Output vector of net production rates. + * Length: m_kk. + */ + virtual void getNetProductionRates(doublereal* wdot) { + err("getNetProductionRates"); + } - //@} + //@} - /** - * @name Reaction Mechanism Informational Query Routines - */ - //@{ + /** + * @name Reaction Mechanism Informational Query Routines + */ + //@{ - /** - * Stoichiometric coefficient of species k as a reactant in - * reaction i. - * - * @param k kinetic species index - * @param i reaction index - */ - virtual doublereal reactantStoichCoeff(int k, int i) const { - err("reactantStoichCoeff"); - return -1.0; - } + /** + * Stoichiometric coefficient of species k as a reactant in + * reaction i. + * + * @param k kinetic species index + * @param i reaction index + */ + virtual doublereal reactantStoichCoeff(int k, int i) const { + err("reactantStoichCoeff"); + return -1.0; + } - /** - * Stoichiometric coefficient of species k as a product in - * reaction i. - * - * @param k kinetic species index - * @param i reaction index - */ - virtual doublereal productStoichCoeff(int k, int i) const { - err("productStoichCoeff"); - return -1.0; - } + /** + * Stoichiometric coefficient of species k as a product in + * reaction i. + * + * @param k kinetic species index + * @param i reaction index + */ + virtual doublereal productStoichCoeff(int k, int i) const { + err("productStoichCoeff"); + return -1.0; + } - /** - * reactant Order of species k in reaction i. - * - * @param k kinetic species index - * @param i reaction index - */ - virtual doublereal reactantOrder(int k, int i) const { - err("reactantOrder"); - return -1.0; - } + /** + * reactant Order of species k in reaction i. + * + * @param k kinetic species index + * @param i reaction index + */ + virtual doublereal reactantOrder(int k, int i) const { + err("reactantOrder"); + return -1.0; + } - /** - * Returns a read-only reference to the vector of reactant - * index numbers for reaction i. - * - * @param i reaction index - */ - virtual const vector_int& reactants(int i) const { - return m_reactants[i]; - } + /** + * Returns a read-only reference to the vector of reactant + * index numbers for reaction i. + * + * @param i reaction index + */ + virtual const vector_int& reactants(int i) const { + return m_reactants[i]; + } - /** - * Returns a read-only reference to the vector of product - * index numbers for reaction i. - * - * @param i reaction index - */ - virtual const vector_int& products(int i) const { - return m_products[i]; - } + /** + * Returns a read-only reference to the vector of product + * index numbers for reaction i. + * + * @param i reaction index + */ + virtual const vector_int& products(int i) const { + return m_products[i]; + } - /** - * Flag specifying the type of reaction. The legal values and - * their meaning are specific to the particular kinetics - * manager. - * - * @param i reaction index - */ - virtual int reactionType(int i) const { - err("reactionType"); - return -1; - } + /** + * Flag specifying the type of reaction. The legal values and + * their meaning are specific to the particular kinetics + * manager. + * + * @param i reaction index + */ + virtual int reactionType(int i) const { + err("reactionType"); + return -1; + } - /** - * True if reaction i has been declared to be reversible. If - * isReversible(i) is false, then the reverse rate of progress - * for reaction i is always zero. - * - * @param i reaction index - */ - virtual bool isReversible(int i){ - err("isReversible"); - return false; - } + /** + * True if reaction i has been declared to be reversible. If + * isReversible(i) is false, then the reverse rate of progress + * for reaction i is always zero. + * + * @param i reaction index + */ + virtual bool isReversible(int i){ + err("isReversible"); + return false; + } - /** - * Return a std::string representing the reaction. - * - * @param i reaction index - */ - virtual std::string reactionString(int i) const { - err("reactionStd::String"); return ""; - } + /** + * Return a std::string representing the reaction. + * + * @param i reaction index + */ + virtual std::string reactionString(int i) const { + err("reactionStd::String"); return ""; + } - /** - * Return the forward rate constants - * - * length is the number of reactions. units depends - * on many issues. @todo DGG: recommend changing name to - * getFwdRateCoefficients. - * - * @param kfwd Output vector containing the foward reaction rate constants. - * Length: m_ii. - */ - virtual void getFwdRateConstants(doublereal *kfwd) { - err("getFwdRateConstants"); - } + /** + * Return the forward rate constants + * + * length is the number of reactions. units depends + * on many issues. @todo DGG: recommend changing name to + * getFwdRateCoefficients. + * + * @param kfwd Output vector containing the foward reaction rate constants. + * Length: m_ii. + */ + virtual void getFwdRateConstants(doublereal *kfwd) { + err("getFwdRateConstants"); + } - /** - * Return the reverse rate constants. - * - * length is the number of reactions. units depends - * on many issues. Note, this routine will return rate constants - * for irreversible reactions if the default for - * doIrreversible is overridden. @todo DGG: recommend changing name to - * getRevRateCoefficients. - * - * @param krev Output vector of reverse rate constants. - * @param doIrreversible boolean indicating whether irreversible reactions - * should be included. - */ - virtual void getRevRateConstants(doublereal *krev, - bool doIrreversible = false) { - err("getFwdRateConstants"); - } + /** + * Return the reverse rate constants. + * + * length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. @todo DGG: recommend changing name to + * getRevRateCoefficients. + * + * @param krev Output vector of reverse rate constants. + * @param doIrreversible boolean indicating whether irreversible reactions + * should be included. + */ + virtual void getRevRateConstants(doublereal *krev, + bool doIrreversible = false) { + err("getFwdRateConstants"); + } - /** - * Return the activation energies in Kelvin. - * - * length is the number of reactions - * - * @param E Ouptut vector of activation energies. - * Length: m_ii. - */ - virtual void getActivationEnergies(doublereal *E) { - err("getActivationEnergies"); - } + /** + * Return the activation energies in Kelvin. + * + * length is the number of reactions + * + * @param E Ouptut vector of activation energies. + * Length: m_ii. + */ + virtual void getActivationEnergies(doublereal *E) { + err("getActivationEnergies"); + } - //@} - /** - * @name Reaction Mechanism Construction - */ - //@{ + //@} + /** + * @name Reaction Mechanism Construction + */ + //@{ - //! Add a phase to the kinetics manager object. - /*! - * This must be done before the function init() is called or - * before any reactions are input. - * The following fields are updated: - * m_start -> vector of integers, containing the - * starting position of the species for - * each phase in the kinetics mechanism. - * m_surfphase -> index of the surface phase. - * m_thermo -> vector of pointers to ThermoPhase phases - * that participate in the kinetics - * mechanism. - * m_phaseindex -> map containing the std::string id of each - * ThermoPhase phase as a key and the - * index of the phase within the kinetics - * manager object as the value. - * - * @param thermo Reference to the ThermoPhase to be added. - */ - virtual void addPhase(thermo_t& thermo); + //! Add a phase to the kinetics manager object. + /*! + * This must be done before the function init() is called or + * before any reactions are input. + * The following fields are updated: + * m_start -> vector of integers, containing the + * starting position of the species for + * each phase in the kinetics mechanism. + * m_surfphase -> index of the surface phase. + * m_thermo -> vector of pointers to ThermoPhase phases + * that participate in the kinetics + * mechanism. + * m_phaseindex -> map containing the std::string id of each + * ThermoPhase phase as a key and the + * index of the phase within the kinetics + * manager object as the value. + * + * @param thermo Reference to the ThermoPhase to be added. + */ + virtual void addPhase(thermo_t& thermo); - /** - * Prepare the class for the addition of reactions. This - * method is called by function importKinetics after all - * phases have been added but before any reactions have - * been. The base class method does nothing, but derived - * classes may use this to perform any initialization - * (allocating arrays, etc.) that requires knowing the phases - * and species, but before any reactions are added. - */ - virtual void init() {} + /** + * Prepare the class for the addition of reactions. This + * method is called by function importKinetics after all + * phases have been added but before any reactions have + * been. The base class method does nothing, but derived + * classes may use this to perform any initialization + * (allocating arrays, etc.) that requires knowing the phases + * and species, but before any reactions are added. + */ + virtual void init() {} - /** - * Finish adding reactions and prepare for use. This method is - * called by function importKinetics after all reactions have - * been entered into the mechanism and before the mechanism is - * used to calculate reaction rates. The base class method - * does nothing, but derived classes may use this to perform - * any initialization (allocating arrays, etc.) that must be - * done after the reactions are entered. - */ - virtual void finalize() {} + /** + * Finish adding reactions and prepare for use. This method is + * called by function importKinetics after all reactions have + * been entered into the mechanism and before the mechanism is + * used to calculate reaction rates. The base class method + * does nothing, but derived classes may use this to perform + * any initialization (allocating arrays, etc.) that must be + * done after the reactions are entered. + */ + virtual void finalize() {} - /** - * Add a single reaction to the mechanism. This routine - * must be called after init() and before finalize(). - * - * @param r Reference to the ReactionData object for the reaction - * to be added. - */ - virtual void addReaction(const ReactionData& r) { - err("addReaction"); - } + /** + * Add a single reaction to the mechanism. This routine + * must be called after init() and before finalize(). + * + * @param r Reference to the ReactionData object for the reaction + * to be added. + */ + virtual void addReaction(const ReactionData& r) { + err("addReaction"); + } - virtual const std::vector& reactantGroups(int i) { - //err("reactantGroups"); - return m_dummygroups; - } + virtual const std::vector& reactantGroups(int i) { + //err("reactantGroups"); + return m_dummygroups; + } - virtual const std::vector& productGroups(int i) { - //err("productGroups"); - return m_dummygroups; - } + virtual const std::vector& productGroups(int i) { + //err("productGroups"); + return m_dummygroups; + } - //@} - /** - * @name Altering Reaction Rates - * - * These methods alter reaction rates. They are designed - * primarily for carrying out sensitivity analysis, but may be - * used for any purpose requiring dynamic alteration of rate - * constants. For each reaction, a real-valued multiplier may - * be defined that multiplies the reaction rate - * coefficient. The multiplier may be set to zero to - * completely remove a reaction from the mechanism. - */ - //@{ + //@} + /** + * @name Altering Reaction Rates + * + * These methods alter reaction rates. They are designed + * primarily for carrying out sensitivity analysis, but may be + * used for any purpose requiring dynamic alteration of rate + * constants. For each reaction, a real-valued multiplier may + * be defined that multiplies the reaction rate + * coefficient. The multiplier may be set to zero to + * completely remove a reaction from the mechanism. + */ + //@{ - /// The current value of the multiplier for reaction i. - /*! - * @param i index of the reaction - */ - doublereal multiplier(int i) const {return m_perturb[i];} + /// The current value of the multiplier for reaction i. + /*! + * @param i index of the reaction + */ + doublereal multiplier(int i) const {return m_perturb[i];} - /// Set the multiplier for reaction i to f. - /*! - * @param i index of the reaction - * @param f value of the multiplier. - */ - void setMultiplier(int i, doublereal f) {m_perturb[i] = f;} + /// Set the multiplier for reaction i to f. + /*! + * @param i index of the reaction + * @param f value of the multiplier. + */ + void setMultiplier(int i, doublereal f) {m_perturb[i] = f;} - //@} + //@} - /** - * Increment the number of reactions in the mechanism by one. - * @todo Should be protected? - */ - void incrementRxnCount() { m_ii++; m_perturb.push_back(1.0); } + /** + * Increment the number of reactions in the mechanism by one. + * @todo Should be protected? + */ + void incrementRxnCount() { m_ii++; m_perturb.push_back(1.0); } - /** - * Returns true if the kinetics manager has been properly - * initialized and finalized. - */ - virtual bool ready() const { - return false; - } + /** + * Returns true if the kinetics manager has been properly + * initialized and finalized. + */ + virtual bool ready() const { + return false; + } - /** - * Extract from array \c data the portion pertaining to phase \c phase. - * - * @param data data - * @param phase phase - * @param phase_data phase_data - */ - void selectPhase(const doublereal* data, const thermo_t* phase, - doublereal* phase_data); + /** + * Extract from array \c data the portion pertaining to phase \c phase. + * + * @param data data + * @param phase phase + * @param phase_data phase_data + */ + void selectPhase(const doublereal* data, const thermo_t* phase, + doublereal* phase_data); - /// For internal use. May be removed in a future release. - int index(){ return m_index; } + /// For internal use. May be removed in a future release. + int index(){ return m_index; } - //! Set the index of the Kinetics Manager - /*! - * @param index input index - */ - void setIndex(int index) { m_index = index; } + //! Set the index of the Kinetics Manager + /*! + * @param index input index + */ + void setIndex(int index) { m_index = index; } - protected: + protected: - //! Number of reactions in the mechanism - int m_ii; + //! Number of reactions in the mechanism + int m_ii; - /// Vector of perturbation factors for each reaction's rate of - /// progress vector. It is initialized to one. - /// - vector_fp m_perturb; + /// Vector of perturbation factors for each reaction's rate of + /// progress vector. It is initialized to one. + /// + vector_fp m_perturb; - /** - * This is a vector of vectors containing the reactants for - * each reaction. The outer vector is over the number of - * reactions, m_ii. The inner vector is a list of species - * indices. If the stoichiometric coefficient for a reactant - * is greater than one, then the reactant is listed - * contiguously in the vector a number of times equal to its - * stoichiometric coefficient. - * NOTE: These vectors will be wrong if there are real - * stoichiometric coefficients in the expression. - */ - std::vector m_reactants; + /** + * This is a vector of vectors containing the reactants for + * each reaction. The outer vector is over the number of + * reactions, m_ii. The inner vector is a list of species + * indices. If the stoichiometric coefficient for a reactant + * is greater than one, then the reactant is listed + * contiguously in the vector a number of times equal to its + * stoichiometric coefficient. + * NOTE: These vectors will be wrong if there are real + * stoichiometric coefficients in the expression. + */ + std::vector m_reactants; - /** - * This is a vector of vectors containing the products for - * each reaction. The outer vector is over the number of - * reactions, m_ii. The inner vector is a list of species - * indeces. If the stoichiometric coefficient for a product is - * greater than one, then the reactant is listed contiguously - * in the vector a number of times equal to its stoichiometric - * coefficient. - * NOTE: These vectors will be wrong if there are real - * stoichiometric coefficients in the expression. - */ - std::vector m_products; + /** + * This is a vector of vectors containing the products for + * each reaction. The outer vector is over the number of + * reactions, m_ii. The inner vector is a list of species + * indeces. If the stoichiometric coefficient for a product is + * greater than one, then the reactant is listed contiguously + * in the vector a number of times equal to its stoichiometric + * coefficient. + * NOTE: These vectors will be wrong if there are real + * stoichiometric coefficients in the expression. + */ + std::vector m_products; - //! m_thermo is a vector of pointers to ThermoPhase - //! objects. - /*! - * For homogeneous kinetics applications, this vector - * will only have one entry. For interfacial reactions, this - * vector will consist of multiple entries; some of them will - * be surface phases, and the other ones will be bulk phases. - * The order that the objects are listed determines the order - * in which the species comprising each phase are listed in - * the source term vector, originating from the reaction - * mechanism. - * - * Note that this kinetics object doesn't own these ThermoPhase - * objects and is not responsible for creating or deleting - * them. - */ - std::vector m_thermo; + //! m_thermo is a vector of pointers to ThermoPhase + //! objects. + /*! + * For homogeneous kinetics applications, this vector + * will only have one entry. For interfacial reactions, this + * vector will consist of multiple entries; some of them will + * be surface phases, and the other ones will be bulk phases. + * The order that the objects are listed determines the order + * in which the species comprising each phase are listed in + * the source term vector, originating from the reaction + * mechanism. + * + * Note that this kinetics object doesn't own these ThermoPhase + * objects and is not responsible for creating or deleting + * them. + */ + std::vector m_thermo; - /** - * m_start is a vector of integers specifying the beginning position - * for the species vector for the n'th phase in the kinetics - * class. - */ - vector_int m_start; + /** + * m_start is a vector of integers specifying the beginning position + * for the species vector for the n'th phase in the kinetics + * class. + */ + vector_int m_start; - /** - * Mapping of the phase id, i.e., the id attribute in the xml - * phase element to the position of the phase within the - * kinetics object. Positions start with the value of 1. The - * member function, phaseIndex() decrements by one before - * returning the index value, so that missing phases return - * -1. - */ - std::map m_phaseindex; - //! Index of the Kinetics Manager - int m_index; + /** + * Mapping of the phase id, i.e., the id attribute in the xml + * phase element to the position of the phase within the + * kinetics object. Positions start with the value of 1. The + * member function, phaseIndex() decrements by one before + * returning the index value, so that missing phases return + * -1. + */ + std::map m_phaseindex; + //! Index of the Kinetics Manager + int m_index; - /** - * Index in the list of phases of the one surface phase. - */ - int m_surfphase; + /** + * Index in the list of phases of the one surface phase. + */ + int m_surfphase; - /** - * Index in the list of phases of the one phase where the reactions - * occur. - */ - int m_rxnphase; + /** + * Index in the list of phases of the one phase where the reactions + * occur. + */ + int m_rxnphase; - /// number of spatial dimensions of lowest-dimensional phase. - int m_mindim; + /// number of spatial dimensions of lowest-dimensional phase. + int m_mindim; - private: + private: - //! Vector of group lists - std::vector m_dummygroups; + //! Vector of group lists + std::vector m_dummygroups; - //! Function for unhandled situations - /*! - * @param m String error message - */ - void err(std::string m) const; + //! Function for unhandled situations + /*! + * @param m String error message + */ + void err(std::string m) const; - }; + }; //! typedef for the kinetics base class - typedef Kinetics kinetics_t; + typedef Kinetics kinetics_t; }