research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Cleaned up Doxygen documentation for class Kinetics and descendants

Browse source
This commit is contained in:
Ray Speth 2013-04-12 23:05:20 +00:00
parent 50617105b5
commit d103ec8d2f
11 changed files with 409 additions and 1733 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

182
include/cantera/kinetics/AqueousKinetics.h
View file

@ -47,11 +47,10 @@ class AqueousKinetics : public Kinetics
public:
/**
* @name Constructors and General Information
*/
//@{
/// Constructor.
//! @name Constructors
//! @{
/// Constructor. Creates an empty reaction mechanism.
AqueousKinetics(thermo_t* thermo = 0);
AqueousKinetics(const AqueousKinetics& right);
@ -74,6 +73,7 @@ public:
* m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
//@}
virtual int type() const {
return cAqueousKinetics;
@ -87,166 +87,56 @@ public:
return m_prxn[k][i];
}
//@}
/**
* @name Reaction Rates Of Progress
*/
//! @name Reaction Rates Of Progress
//@{
/**
* 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.
*/
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
updateROP();
std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
}
/**
* Reverse rates of progress.
* Return the reverse rates of progress in array revROP, which
* must be dimensioned at least as large as the total number
* of reactions.
*/
virtual void getRevRatesOfProgress(doublereal* revROP) {
updateROP();
std::copy(m_ropr.begin(), m_ropr.end(), revROP);
}
/**
* 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.
*/
virtual void getNetRatesOfProgress(doublereal* netROP) {
updateROP();
std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
}
/**
* 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.
*/
virtual void getEquilibriumConstants(doublereal* kc);
/**
* Return the array of values for the reaction gibbs free energy
* change.
* These values depend on the species concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaGibbs(doublereal* deltaG);
/**
* Return the array of values for the reaction enthalpy change.
* These values depend upon the species concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaEnthalpy(doublereal* deltaH);
/**
* Return the array of values for the reactions change in
* entropy.
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaEntropy(doublereal* deltaS);
/**
* Return the array of values for the reaction
* standard state Gibbs free energy change.
* These values do not depend on the species
* concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaSSGibbs(doublereal* deltaG);
/**
* Return the array of values for the change in the
* standard state enthalpies of reaction.
* These values do not depend upon the concentration
* of the solution.
*
* units = J kmol-1
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
/**
* Return the array of values for the change in the
* standard state entropies for each reaction.
* These values do not depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaSSEntropy(doublereal* deltaS);
//@}
/**
* @name Species Production Rates
*/
//@{
//! @}
//! @name Species Production Rates
//! @{
//! Return the species net production rates
/*!
* Species net production rates [kmol/m^3/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.
*
* @param net Array of species production rates.
* units kmol m-3 s-1
*/
virtual void getNetProductionRates(doublereal* net) {
updateROP();
m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net);
}
/**
* Species creation rates [kmol/m^3]. Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species.
*
*/
virtual void getCreationRates(doublereal* cdot) {
updateROP();
m_rxnstoich.getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot);
}
/**
* Species destruction rates [kmol/m^3]. Return the species
* destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species.
*
*/
virtual void getDestructionRates(doublereal* ddot) {
updateROP();
m_rxnstoich.getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot);
}
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
//! @}
//! @name Reaction Mechanism Informational Query Routines
//! @{
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(size_t i) const {
return m_index[i].first;
}
@ -255,11 +145,6 @@ public:
return m_rxneqn[i];
}
/**
* 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.
*/
virtual bool isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
@ -269,45 +154,25 @@ public:
}
}
/**
* Return the forward rate constants
*
* length is the number of reactions. units depends
* on many issues.
*/
virtual void getFwdRateConstants(doublereal* kfwd);
/**
* 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.
*/
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
//@}
/**
* @name Reaction Mechanism Setup Routines
*/
//@{
//! @}
//! @name Reaction Mechanism Setup Routines
//! @{
virtual void init();
/// Add a reaction to the mechanism.
virtual void addReaction(ReactionData& r);
virtual void finalize();
virtual bool ready() const;
virtual void update_T();
virtual void update_C();
void updateROP();
const std::vector<grouplist_t>& reactantGroups(size_t i) {
return m_rgroups[i];
}
@ -315,8 +180,16 @@ public:
return m_pgroups[i];
}
/*!
* Update temperature-dependent portions of reaction rates and
* falloff functions.
*/
void _update_rates_T();
/*!
* Update properties that depend on concentrations. Currently only
* the enhanced collision partner concentrations are updated here.
*/
void _update_rates_C();
//@}
@ -382,11 +255,14 @@ private:
void addElementaryReaction(ReactionData& r);
void installReagents(const ReactionData& r);
void installGroups(size_t irxn, const std::vector<grouplist_t>& r,
const std::vector<grouplist_t>& p);
/**
* Update the equilibrium constants in molar units.
*/
void updateKc();
void registerReaction(size_t rxnNumber, int type, size_t loc) {

24
include/cantera/kinetics/EdgeKinetics.h
View file

@ -6,7 +6,6 @@
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_EDGEKINETICS_H
#define CT_EDGEKINETICS_H
@ -14,20 +13,14 @@
namespace Cantera
{
/**
* Heterogeneous reactions at one-dimensional interfaces between
* multiple adjacent two-dimensional surfaces.
*/
class EdgeKinetics : public InterfaceKinetics
{
public:
/**
* Constructor
*
*/
//! Constructor
EdgeKinetics() : InterfaceKinetics() {}
/// Destructor.
@ -45,27 +38,12 @@ public:
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.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const {
EdgeKinetics* iK = new EdgeKinetics(*this);
iK->assignShallowPointers(tpVector);
return iK;
}
/**
* Identifies the subclass of the Kinetics manager type.
* These are listed in mix_defs.h.
*/
virtual int type() const {
return cEdgeKinetics;
}

10
include/cantera/kinetics/GRI_30_Kinetics.h
View file

@ -1,5 +1,4 @@
/**
*
* @file GRI_30_Kinetics.h
*
* @ingroup chemkinetics
@ -7,7 +6,6 @@
// Copyright 2001 California Institute of Technology
#ifndef CT_GRI30_KINETICS_H
#define CT_GRI30_KINETICS_H
@ -15,7 +13,6 @@
namespace Cantera
{
const int cGRI_30_Kinetics = cGasKinetics + 1;
/**
@ -23,9 +20,7 @@ const int cGRI_30_Kinetics = cGasKinetics + 1;
*/
class GRI_30_Kinetics : public GasKinetics
{
public:
/// Default constructor.
GRI_30_Kinetics(thermo_t* th=0);
@ -44,7 +39,12 @@ public:
private:
void gri30_update_rates_T();
void gri30_updateROP();
/**
* Update the equilibrium constants in molar units.
*/
void gri30_updateKc();
void get_wdot(const doublereal* rop, doublereal* wdot);
void update_kc(const doublereal* grt, doublereal c0, doublereal* rkc);
void update_rates(doublereal t, doublereal tlog, doublereal* rf);

219
include/cantera/kinetics/GasKinetics.h
View file

@ -6,7 +6,6 @@
// Copyright 2001 California Institute of Technology
#ifndef CT_GASKINETICS_H
#define CT_GASKINETICS_H
@ -41,13 +40,9 @@ class ReactionData;
*/
class GasKinetics : public Kinetics
{
public:
/**
* @name Constructors and General Information
*/
//@{
//! @name Constructors and General Information
//! @{
//! Constructor.
/*!
@ -55,47 +50,17 @@ public:
*/
GasKinetics(thermo_t* thermo = 0);
//!Copy Constructor for the %GasKinetics object.
/*!
* Currently, this is not fully implemented. If called it will
* throw an exception.
*
* @param right object to be copied
*/
//! Copy Constructor
GasKinetics(const GasKinetics& right);
//! Destructor.
virtual ~GasKinetics();
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %GasKinetics object to be copied into the
* current one.
*/
GasKinetics& operator=(const GasKinetics& right);
//! 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.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
//! 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 {
return cGasKinetics;
}
@ -108,165 +73,46 @@ public:
return m_prxn[k][i];
}
//@}
/**
* @name Reaction Rates Of Progress
*/
//@{
/**
* 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.
*/
//! @}
//! @name Reaction Rates Of Progress
//! @{
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
updateROP();
std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
}
/**
* Reverse rates of progress.
* Return the reverse rates of progress in array revROP, which
* must be dimensioned at least as large as the total number
* of reactions.
*/
virtual void getRevRatesOfProgress(doublereal* revROP) {
updateROP();
std::copy(m_ropr.begin(), m_ropr.end(), revROP);
}
/**
* 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.
*/
virtual void getNetRatesOfProgress(doublereal* netROP) {
updateROP();
std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
}
/**
* 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.
*/
virtual void getEquilibriumConstants(doublereal* kc);
/**
* Return the array of values for the reaction gibbs free energy
* change.
* These values depend on the species concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaGibbs(doublereal* deltaG);
/**
* Return the array of values for the reaction enthalpy change.
* These values depend upon the species concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaEnthalpy(doublereal* deltaH);
/**
* Return the array of values for the reactions change in
* entropy.
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaEntropy(doublereal* deltaS);
/**
* Return the array of values for the reaction
* standard state Gibbs free energy change.
* These values do not depend on the species
* concentrations.
*
* units = J kmol-1
*/
virtual void getDeltaSSGibbs(doublereal* deltaG);
/**
* Return the array of values for the change in the
* standard state enthalpies of reaction.
* These values do not depend upon the concentration
* of the solution.
*
* units = J kmol-1
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
/**
* Return the array of values for the change in the
* standard state entropies for each reaction.
* These values do not depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaSSEntropy(doublereal* deltaS);
//@}
/**
* @name Species Production Rates
*/
//@{
//! @}
//! @name Species Production Rates
//! @{
//! Return the species net production rates
/*!
* Species net production rates [kmol/m^3/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.
*
* @param net Array of species production rates.
* units kmol m-3 s-1
*/
virtual void getNetProductionRates(doublereal* net);
//! Return the species creation rates
/*!
* Species creation rates [kmol/m^3]. Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species.
*
* @param cdot Array of species creation rates.
* units kmol m-3 s-1
*/
virtual void getCreationRates(doublereal* cdot);
//! Return a vector of the species destruction rates
/*!
* Species destruction rates [kmol/m^3]. Return the species
* destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species.
*
*
* @param ddot Array of species destruction rates.
* units kmol m-3 s-1
*
*/
virtual void getDestructionRates(doublereal* ddot);
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
//! @}
//! @name Reaction Mechanism Informational Query Routines
//! @{
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(size_t i) const {
return m_index[i].first;
}
@ -275,11 +121,6 @@ public:
return m_rxneqn[i];
}
/**
* 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.
*/
virtual bool isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
@ -289,38 +130,19 @@ public:
}
}
/**
* Return the forward rate constants
*
* length is the number of reactions. units depends
* on many issues.
*/
virtual void getFwdRateConstants(doublereal* kfwd);
/**
* 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.
*/
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
//@}
/**
* @name Reaction Mechanism Setup Routines
*/
//@{
//! @}
//! @name Reaction Mechanism Setup Routines
//! @{
virtual void init();
/// Add a reaction to the mechanism.
virtual void addReaction(ReactionData& r);
virtual void finalize();
virtual bool ready() const;
//@}
void updateROP();
@ -341,10 +163,7 @@ public:
//! reactions.
virtual void update_rates_C();
//@}
protected:
size_t m_nfall;
std::vector<size_t> m_fallindx;
@ -415,9 +234,7 @@ protected:
void processFalloffReactions();
vector_fp m_grt;
private:
size_t reactionNumber() {
return m_ii;
}
@ -433,6 +250,8 @@ private:
void installGroups(size_t irxn, const std::vector<grouplist_t>& r,
const std::vector<grouplist_t>& p);
//! Update the equilibrium constants in molar units.
void updateKc();
void registerReaction(size_t rxnNumber, int type_, size_t loc) {

540
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -33,40 +33,42 @@ class ImplicitSurfChem;
//! A kinetics manager for heterogeneous reaction mechanisms. The
//! reactions are assumed to occur at a 2D interface between two 3D phases.
/*!
* There are some important additions to the behavior of the kinetics class
* due to the presence of multiple phases and a heterogeneous interface. If
* a reactant phase doesn't exists, i.e., has a mole number of zero, a
* heterogeneous reaction can not proceed from reactants to products. Note it
* could perhaps proceed from products to reactants if all of the product
* phases exist.
*
* There are some important additions to the behavior of the kinetics class due to the
* presence of multiple phases and a heterogeneous interface. If a reactant phase
* doesn't exists, i.e., has a mole number of zero, a heterogeneous reaction can not
* proceed from reactants to products. Note it could perhaps proceed from products to
* reactants if all of the product phases exist.
* In order to make the determination of whether a phase exists or not
* actually involves the specification of additional information to the
* kinetics object., which heretofore has only had access to intrinsic field
* information about the phases (i.e., temperature pressure, and mole
* fraction).
*
* In order to make the determination of whether a phase exists or not actually involves
* the specification of additional information to the kinetics object., which heretofore
* has only had access to intrinsic field information about the phases (i.e., temperature
* pressure, and mole fraction).
* The extrinsic specification of whether a phase exists or not must be
* specified on top of the intrinsic calculation of the reaction rate. This
* class carries a set of booleans indicating whether a phase in the
* heterogeneous mechanism exists or not.
*
* The extrinsic specification of whether a phase exists or not must be specified on top of the
* intrinsic calculation of the reaction rate. This routine carries a set of
* booleans indicating whether a phase in the heterogeneous mechanism exists or not.
*
* Additionally, the routine carries a set of booleans around indicating whether a product
* phase is stable or not. If a phase is not thermodynamically stable, it may be the case that
* a particular reaction in a heterogeneous mechanism will create a product species in the
* unstable phase. However, other reactions in the mechanism will destruct that species.
* This may cause oscillations in the formation of the unstable phase from time step to time
* step within a ODE solver, in practice. In order to avoid this situation, a set of
* booleans is tracked which sets the stability of a phase. If a phase is deemed to be unstable,
* then species in that phase will not be allowed to be birthed by the kinetics operator.
* Nonexistent phases are deemed to be unstable by default, but this can be changed.
* Additionally, the class carries a set of booleans around indicating
* whether a product phase is stable or not. If a phase is not
* thermodynamically stable, it may be the case that a particular reaction in
* a heterogeneous mechanism will create a product species in the unstable
* phase. However, other reactions in the mechanism will destruct that
* species. This may cause oscillations in the formation of the unstable
* phase from time step to time step within a ODE solver, in practice. In
* order to avoid this situation, a set of booleans is tracked which sets the
* stability of a phase. If a phase is deemed to be unstable, then species in
* that phase will not be allowed to be birthed by the kinetics operator.
* Nonexistent phases are deemed to be unstable by default, but this can be
* changed.
*
* @ingroup chemkinetics
*/
class InterfaceKinetics : public Kinetics
{
public:
//! Constructor
/*!
* @param thermo The optional parameter may be used to initialize
@ -78,42 +80,17 @@ public:
*/
InterfaceKinetics(thermo_t* thermo = 0);
/// Destructor.
virtual ~InterfaceKinetics();
//! 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 %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.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
//! Return the type of the kinetics object
virtual int type() const;
//! Set the electric potential in the nth phase
@ -123,204 +100,51 @@ public:
*/
void setElectricPotential(int n, doublereal V);
///
/// @name Reaction Rates Of Progress
///
//@{
//! @name Reaction Rates Of Progress
//! @{
//! Return the forward rates of progress for each reaction
/*!
* @param fwdROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getFwdRatesOfProgress(doublereal* fwdROP);
//! Return the reverse rates of progress for each reaction
/*!
* @param revROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getRevRatesOfProgress(doublereal* revROP);
//! Return the net rates of progress for each reaction
/*!
* @param netROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getNetRatesOfProgress(doublereal* netROP);
//! Get the equilibrium constants of all reactions, whether
//! the reaction is reversible or not.
/*!
* @param kc Returns the concentration equation constant for the reaction.
* Length is the number of reactions
*/
virtual void getEquilibriumConstants(doublereal* kc);
void getExchangeCurrentQuantities();
//! 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);
//! Return the vector of values for the reaction electrochemical free energy change.
/*!
* These values depend upon the concentration of the solution and
* the voltage of the phases
*
* units = J kmol-1
*
* @param deltaM Output vector of deltaM's for reactions
* Length: m_ii.
*/
virtual void getDeltaElectrochemPotentials(doublereal* deltaM);
/**
* Return the vector of values for the reactions change in
* enthalpy.
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1
*/
virtual void getDeltaEnthalpy(doublereal* deltaH);
//! Return the vector of values for the change in
//! entropy due to each reaction
/*!
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*
* @param deltaS vector of Enthalpy changes
* Length = m_ii, number of reactions
*
*/
virtual void getDeltaEntropy(doublereal* deltaS);
//! 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.
*
* @param deltaG vector of rxn SS free energy changes
* units = J kmol-1
*/
virtual void getDeltaSSGibbs(doublereal* deltaG);
//! 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.
*
* @param deltaH vector of rxn SS enthalpy changes
* units = J kmol-1
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
//! 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.
*
* @param deltaS vector of rxn SS entropy changes
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaSSEntropy(doublereal* deltaS);
//! @}
//! @name Species Production Rates
//! @{
//@}
/**
* @name Species Production Rates
*/
//@{
//! Returns the Species creation rates [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 of the kinetics
* model
*
* @param cdot Vector containing creation rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getCreationRates(doublereal* cdot);
//! Return the Species destruction rates [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 in all phases of the kinetics model
*
* @param ddot Vector containing destruction rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getDestructionRates(doublereal* ddot);
//! Return the species net production rates [kmol/m^2/s].
/*!
* Species net production rates [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 in all phases of the kinetics
* model
*
* @param net Vector of species production rates.
* units kmol m-d s-1, where d is dimension.
*/
virtual void getNetProductionRates(doublereal* net);
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
//! @}
//! @name Reaction Mechanism Informational Query Routines
//! @{
/**
* Stoichiometric coefficient of species k as a reactant in
* reaction i.
*/
virtual doublereal reactantStoichCoeff(size_t k, size_t i) const {
return m_rrxn[k][i];
}
/**
* Stoichiometric coefficient of species k as a product in
* reaction i.
*/
virtual doublereal productStoichCoeff(size_t k, size_t i) const {
return m_prxn[k][i];
}
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(size_t i) const {
return m_index[i].first;
}
//! Get the vector of activity concentrations used in the kinetics object
/*!
* @param conc (output) Vector of activity concentrations. Length is
* equal to the number of species in the kinetics object
*/
virtual void getActivityConcentrations(doublereal* const conc);
//! Return the charge transfer rxn Beta parameter for the ith reaction
@ -332,20 +156,13 @@ public:
* no information is known, as a value of 0.5 pertains to a
* symmetric transition state. The value can vary between 0 to 1.
*
*
* @param irxn Reaction number in the kinetics mechanism
*
* @return
* Beta parameter. This defaults to zero, even for charge transfer
* reactions.
* @return Beta parameter. This defaults to zero, even for charge
* transfer reactions.
*/
doublereal electrochem_beta(size_t irxn) const;
/**
* 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.
*/
virtual bool isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
@ -355,34 +172,26 @@ public:
}
}
/**
* Return a string representing the reaction.
*/
virtual std::string reactionString(size_t i) const {
return m_rxneqn[i];
}
virtual void getFwdRateConstants(doublereal* kfwd);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
virtual void getActivationEnergies(doublereal* E);
//@}
/**
* @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.
*
* This function calls the Kinetics operator addPhase.
* It also sets the following functions
* This function calls Kinetics::addPhase(). It also sets the following
* fields:
*
* m_phaseExists[]
*
@ -390,33 +199,11 @@ public:
*/
virtual void addPhase(thermo_t& thermo);
//! Prepare the class for the addition of reactions.
/*!
* This function must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
virtual void init();
//! Add a single reaction to the mechanism.
/*!
* @param r Reference to a ReactionData object containing all of
* the info needed to describe the reaction.
*/
virtual void addReaction(ReactionData& r);
//! Finish adding reactions and prepare for use.
/*!
* This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*/
virtual void finalize();
virtual bool ready() const;
//! @}
//! Internal routine that updates the Rates of Progress of the reactions
/*!
@ -424,13 +211,9 @@ public:
*/
void updateROP();
//! Update properties that depend on temperature
/*!
* This is called to update all of the properties that depend on temperature
*
* Current objects that this function updates
* Current objects that this function updates:
* m_kdata->m_logtemp
* m_kdata->m_rfn
* m_rates.
@ -439,14 +222,15 @@ public:
void _update_rates_T();
//! Update properties that depend on the electric potential
/*!
* This is called to update all of the properties that depend on potential
*/
void _update_rates_phi();
//! Update properties that depend on the species mole fractions and/or concentration
//! Update properties that depend on the species mole fractions and/or
//! concentration,
/*!
* This is called to update all of the properties that depend on concentration
* This method fills out the array of generalized concentrations by
* calling method getActivityConcentrations for each phase, which classes
* representing phases should overload to return the appropriate
* quantities.
*/
void _update_rates_C();
@ -459,14 +243,12 @@ public:
* \dot {\theta}_k = \dot s_k (\sigma_k / s_0)
* \f]
*
*
* @param tstep Time value to advance the surface coverages
*/
void advanceCoverages(doublereal tstep);
//! Solve for the pseudo steady-state of the surface problem
/*!
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
@ -497,7 +279,6 @@ public:
void checkPartialEquil();
size_t reactionNumber() const {
return m_ii;
}
@ -506,6 +287,12 @@ public:
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
/**
* Update the equilibrium constants in molar units for all reversible
* reactions. Irreversible reactions have their equilibrium constant set
* to zero. For reactions involving charged species the equilibrium
* constant is adjusted according to the electrostatic potential.
*/
void updateKc();
//! Write values into m_index
@ -531,36 +318,37 @@ public:
//! When an electrode reaction rate is optionally specified in terms of its
//! exchange current density, extra vectors need to be precalculated
/*!
*
*/
void applyExchangeCurrentDensityFormulation(doublereal* const kfwd);
//! Set the existence of a phase in the reaction object
/*!
* Tell the kinetics object whether a phase in the object exists.
* This is actually an extrinsic specification that must be carried out on top of the
* intrinsic calculation of the reaction rate.
* The routine will also flip the IsStable boolean within the kinetics object as well.
* Tell the kinetics object whether a phase in the object exists. This is
* actually an extrinsic specification that must be carried out on top of
* the intrinsic calculation of the reaction rate. The routine will also
* flip the IsStable boolean within the kinetics object as well.
*
* @param iphase Index of the phase. This is the order within the internal thermo vector object
* @param iphase Index of the phase. This is the order within the
* internal thermo vector object
* @param exists Boolean indicating whether the phase exists or not
*/
void setPhaseExistence(const size_t iphase, const int exists);
//! Set the stability of a phase in the reaction object
/*!
* Tell the kinetics object whether a phase in the object is stable. Species in an unstable phase
* will not be allowed to have a positive rate of formation from this kinetics object.
* This is actually an extrinsic specification that must be carried out on top of the
* intrinsic calculation of the reaction rate.
* Tell the kinetics object whether a phase in the object is stable.
* Species in an unstable phase will not be allowed to have a positive
* rate of formation from this kinetics object. This is actually an
* extrinsic specification that must be carried out on top of the
* intrinsic calculation of the reaction rate.
*
* While conceptually not needed since kinetics is consistent with thermo when taken as a whole,
* in practice it has found to be very useful to turn off the creation of phases which shouldn't
* be forming. Typically this can reduce the oscillations in phase formation and destruction
* which are observed.
* While conceptually not needed since kinetics is consistent with thermo
* when taken as a whole, in practice it has found to be very useful to
* turn off the creation of phases which shouldn't be forming. Typically
* this can reduce the oscillations in phase formation and destruction
* which are observed.
*
* @param iphase Index of the phase. This is the order within the internal thermo vector object
* @param iphase Index of the phase. This is the order within the
* internal thermo vector object
* @param isStable Flag indicating whether the phase is stable or not
*/
void setPhaseStability(const size_t iphase, const int isStable);
@ -568,34 +356,30 @@ public:
//! Gets the phase existence int for the ith phase
/*!
* @param iphase Phase Id
*
* @return Returns the int specifying whether the kinetics object thinks the phase exists
* or not. If it exists, then species in that phase can be a reactant in reactions.
* @return The int specifying whether the kinetics object thinks the phase
* exists or not. If it exists, then species in that phase can be
* a reactant in reactions.
*/
int phaseExistence(const size_t iphase) const;
//! Gets the phase stability int for the ith phase
/*!
* @param iphase Phase Id
*
* @return Returns the int specifying whether the kinetics object thinks the phase is stable
* with nonzero mole numbers.
* If it stable, then the kinetics object will allow for rates of production of
* of species in that phase that are positive.
* @return The int specifying whether the kinetics object thinks the phase
* is stable with nonzero mole numbers. If it stable, then the
* kinetics object will allow for rates of production of of
* species in that phase that are positive.
*/
int phaseStability(const size_t iphase) const;
protected:
//! Temporary work vector of length m_kk
vector_fp m_grt;
//! List of reactions numbers which are reversible reactions
/*!
* This is a vector of reaction numbers. Each reaction
* in the list is reversible.
* Length = number of reversible reactions
* This is a vector of reaction numbers. Each reaction in the list is
* reversible. Length = number of reversible reactions
*/
std::vector<size_t> m_revindex;
@ -606,11 +390,10 @@ protected:
*/
Rate1<SurfaceArrhenius> m_rates;
bool m_redo_rates;
bool m_redo_rates;
/**
* Vector of information about reactions in the
* mechanism.
* Vector of information about reactions in the mechanism.
* The key is the reaction index (0 < i < m_ii).
* The first pair is the reactionType of the reaction.
* The second pair is ...
@ -619,17 +402,15 @@ protected:
//! Vector of irreversible reaction numbers
/*!
* vector containing the reaction numbers of irreversible
* reactions.
* vector containing the reaction numbers of irreversible reactions.
*/
std::vector<size_t> m_irrev;
//! Stoichiometric manager for the reaction mechanism
/*!
* This is the manager for the kinetics mechanism that
* handles turning reaction extents into species
* production rates and also handles turning thermo
* properties into reaction thermo properties.
* This is the manager for the kinetics mechanism that handles turning
* reaction extents into species production rates and also handles
* turning thermo properties into reaction thermo properties.
*/
ReactionStoichMgr m_rxnstoich;
@ -639,15 +420,14 @@ protected:
//! Number of reversible reactions in the mechanism
size_t m_nrev;
//! m_rrxn is a vector of maps, containing the reactant
//! stoichiometric coefficient information
/*!
* m_rrxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* reactant stoichiometric coefficient for the species being the value.
* m_rrxn has a length equal to the total number of species in the
* kinetics object. For each species, there exists a map, with the
* reaction number being the key, and the reactant stoichiometric
* coefficient for the species being the value.
*
* HKM -> mutable because search sometimes creates extra
* entries. To be fixed in future...
*/
@ -656,57 +436,50 @@ protected:
//! m_prxn is a vector of maps, containing the reactant
//! stoichiometric coefficient information
/**
* m_prxn is a vector of maps. m_prxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* product stoichiometric coefficient for the species being the value.
* m_prxn is a vector of maps. m_prxn has a length equal to the total
* number of species in the kinetics object. For each species, there
* exists a map, with the reaction number being the key, and the product
* stoichiometric coefficient for the species being the value.
*/
mutable std::vector<std::map<size_t, doublereal> > m_prxn;
//! String expression for each rxn
/*!
* Vector of strings of length m_ii, the number of
* reactions, containing the
* string expressions for each reaction
* reactions, containing the string expressions for each reaction
* (e.g., reactants <=> product1 + product2)
*/
std::vector<std::string> m_rxneqn;
//! an array of generalized concentrations for each species
/*!
* An array of generalized concentrations
* \f$ C_k \f$ that are defined such that \f$ a_k = C_k /
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/
* These generalized concentrations are used
* by this kinetics manager class to compute the forward and
* reverse rates of elementary reactions. The "units" for the
* concentrations of each phase depend upon the implementation
* of kinetics within that phase.
* The order of the species within the vector is based on
* the order of listed ThermoPhase objects in the class, and the
* order of the species within each ThermoPhase class.
* An array of generalized concentrations \f$ C_k \f$ that are defined
* such that \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ is a standard
* concentration/ These generalized concentrations are used by this
* kinetics manager class to compute the forward and reverse rates of
* elementary reactions. The "units" for the concentrations of each phase
* depend upon the implementation of kinetics within that phase. The order
* of the species within the vector is based on the order of listed
* ThermoPhase objects in the class, and the order of the species within
* each ThermoPhase class.
*/
vector_fp m_conc;
//! Vector of standard state chemical potentials
/*!
* This vector contains a temporary vector of
* standard state chemical potentials
* for all of the species in the kinetics object
* This vector contains a temporary vector of standard state chemical
* potentials for all of the species in the kinetics object
*
* Length = m_k
* units = J/kmol
* Length = m_k. Units = J/kmol.
*/
vector_fp m_mu0;
//! Vector of phase electric potentials
/*!
* Temporary vector containing the potential of each phase
* in the kinetics object
* Temporary vector containing the potential of each phase in the kinetics
* object.
*
* length = number of phases
* units = Volts
* length = number of phases. Units = Volts.
*/
vector_fp m_phi;
@ -736,17 +509,16 @@ protected:
//! Pointer to the Implicit surface chemistry object
/*!
* Note this object is owned by this InterfaceKinetics
* object. It may only be used to solve this single
* InterfaceKinetics objects's surface problem uncoupled
* from other surface phases.
* Note this object is owned by this InterfaceKinetics object. It may only
* be used to solve this single InterfaceKinetics objects's surface
* problem uncoupled from other surface phases.
*/
ImplicitSurfChem* m_integrator;
vector_fp m_beta;
//! Vector of reaction indexes specifying the id of the current transfer reactions
//! in the mechanism
//! Vector of reaction indexes specifying the id of the current transfer
//! reactions in the mechanism
/*!
* Vector of reaction indices which involve current transfers. This provides
* an index into the m_beta array.
@ -755,8 +527,8 @@ protected:
*/
std::vector<size_t> m_ctrxn;
//! Vector of booleans indicating whether the charge transfer reaction may be
//! described by an exchange current density expression
//! Vector of booleans indicating whether the charge transfer reaction may
//! be described by an exchange current density expression
vector_int m_ctrxn_ecdf;
vector_fp m_StandardConc;
@ -802,80 +574,78 @@ protected:
//! Boolean flag indicating whether any reaction in the mechanism
//! is described by an exchange current density expression
/*!
* If this is true, the standard state gibbs free energy of the reaction and
* the product of the reactant standard concentrations must be precalculated
* in order to calculate the rate constant.
* If this is true, the standard state gibbs free energy of the reaction
* and the product of the reactant standard concentrations must be
* precalculated in order to calculate the rate constant.
*/
bool m_has_exchange_current_density_formulation;
//! Int flag to indicate that some phases in the kinetics mechanism are
//! non-existent.
/*!
* We change the ROP vectors to make sure that non-existent phases are treated
* correctly in the kinetics operator. The value of this is equal to the number
* of phases which don't exist.
* We change the ROP vectors to make sure that non-existent phases are
* treated correctly in the kinetics operator. The value of this is equal
* to the number of phases which don't exist.
*/
int m_phaseExistsCheck;
//! Vector of booleans indicating whether phases exist or not
/*!
* Vector of booleans indicating whether a phase exists or not.
* We use this to set the ROP's so that unphysical things don't happen
* Vector of booleans indicating whether a phase exists or not. We use
* this to set the ROP's so that unphysical things don't happen
*
* length = number of phases in the object
* By default all phases exist.
* length = number of phases in the object. By default all phases exist.
*/
std::vector<bool> m_phaseExists;
//! Vector of int indicating whether phases are stable or not
//! Vector of int indicating whether phases are stable or not
/*!
* Vector of booleans indicating whether a phase is stable or not
* under the current conditions.
* We use this to set the ROP's so that unphysical things don't happen
* Vector of booleans indicating whether a phase is stable or not under
* the current conditions. We use this to set the ROP's so that
* unphysical things don't happen
*
* length = number of phases in the object
* By default all phases are stable
* length = number of phases in the object. By default all phases are stable.
*/
std::vector<int> m_phaseIsStable;
//! Vector of vector of booleans indicating whether a phase participates in a
//! reaction as a reactant
//! Vector of vector of booleans indicating whether a phase participates in a
//! reaction as a reactant
/*!
* m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
* participates in reaction j as a reactant.
* m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
* participates in reaction j as a reactant.
*/
std::vector<std::vector<bool> > m_rxnPhaseIsReactant;
//! Vector of vector of booleans indicating whether a phase participates in a
//! reaction as a product
//! Vector of vector of booleans indicating whether a phase participates in a
//! reaction as a product
/*!
* m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
* participates in reaction j as a product.
* m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
* participates in reaction j as a product.
*/
std::vector<std::vector<bool> > m_rxnPhaseIsProduct;
#ifdef KINETICS_WITH_INTERMEDIATE_ZEROED_PHASES
//! Vector of ints indicating whether zeroed phase is an intermediate for
//! the formation of another phase
//! Vector of ints indicating whether zeroed phase is an intermediate for
//! the formation of another phase
/*!
* If a phase is zeroed out but it is an intermediate, then the phase
* can be formed whether it is stable or not, but the destruction rate of
* species in that phase can't exceed the formation rate for species in that
* phase.
* If a phase is zeroed out but it is an intermediate, then the phase can
* be formed whether it is stable or not, but the destruction rate of
* species in that phase can't exceed the formation rate for species in
* that phase.
*
* length = number of phases in the object
* By default all phases are not intermediates
* length = number of phases in the object. By default all phases are not
* intermediates
*/
std::vector<int> m_phaseIsIntermediate;
int m_numIntermediatePhases;
//! Reaction rate reduction factor for intermediates
//! Reaction rate reduction factor for intermediates
/*!
* Individual reaction rates are reduced to accommodate the requirements of intermediate
* zero phases.
*
* length = number of reactions in the object
* By default all phases are not intermediates
* length = number of reactions in the object. By default all phases are
* not intermediates.
*/
std::vector<doublereal> m_rxnRateFactorPhaseIntermediates;
@ -885,8 +655,6 @@ protected:
#endif
int m_ioFlag;
private:
};
}

377
include/cantera/kinetics/Kinetics.h
View file

@ -143,50 +143,44 @@ public:
/// Destructor.
virtual ~Kinetics();
//!Copy Constructor for the %Kinetics object.
/*!
* Currently, this is not fully implemented. If called it will
* throw an exception.
*/
//!Copy Constructor for the Kinetics object.
Kinetics(const Kinetics&);
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %Kinetics object to be copied into the
* @param right Reference to Kinetics object to be copied into the
* current one.
*/
Kinetics& operator=(const Kinetics& right);
//! 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.
* This function can be used to duplicate objects derived 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.
* These routines are basically wrappers around the derived copy
* constructor.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
* @param tpVector Vector of pointers to ThermoPhase objects. this is the
* #m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
//! Reassign the shallow pointers within the %FKinetics object
//! Reassign the pointers within the Kinetics object
/*!
* This type or routine is absolute necessary because the Kinetics object doesn't
* own the ThermoPhase objects. After a duplication, we need to point to different
* ThermoPhase objects.
* This type or routine is necessary because the Kinetics object doesn't
* own the ThermoPhase objects. After a duplication, we need to point to
* different ThermoPhase objects.
*
* We check that the ThermoPhase objects are aligned in the same order and have
* the following identical properties to the ones that they are replacing.
* id()
* eosType()
* nSpecies()
* the following identical properties to the ones that they are replacing:
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
* - ThermoPhase::id()
* - ThermoPhase::eosType()
* - ThermoPhase::nSpecies()
*
* @param tpVector Vector of pointers to ThermoPhase objects. this is the
* #m_thermo vector within this object
*/
virtual void assignShallowPointers(const std::vector<thermo_t*> & tpVector);
@ -220,12 +214,9 @@ public:
//! Throws an exception if kk is less than nSpecies(). Used before calls
//! which take an array pointer.
void checkSpeciesArraySize(size_t mm) const;
//@}
/**
* @name Information/Lookup Functions about Phases and Species
*/
//! @name Information/Lookup Functions about Phases and Species
//@{
/**
@ -265,36 +256,32 @@ public:
}
/**
* 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 returns -1.
* 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 returns -1.
*/
size_t 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.
* 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.
*/
size_t 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.
* 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.
*/
@ -306,10 +293,9 @@ public:
}
/**
* 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.
* 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.
*/
size_t nTotalSpecies() const {
size_t n=0, np;
@ -346,15 +332,12 @@ public:
return m_start[n] + k;
}
//! Return the std::string name of the kth species in the kinetics
//! manager.
//! Return the 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
* "<unknown>" is returned.
* 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 "<unknown>" is returned.
*
* @param k species index
*/
@ -389,10 +372,9 @@ public:
const 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.
* This function looks up the 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 doesn't match.
*
* @param nm String containing the name of the species.
*/
@ -410,23 +392,17 @@ public:
}
/**
* 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.
* 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
*/
size_t speciesPhaseIndex(size_t k);
//@}
/**
* @name Reaction Rates Of Progress
*/
//@{
//! @}
//! @name Reaction Rates Of Progress
//! @{
//! Return the forward rates of progress of the reactions
/*!
@ -443,9 +419,8 @@ public:
//! 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.
* 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.
@ -455,10 +430,9 @@ public:
}
/**
* 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.
* 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.
*/
@ -466,14 +440,11 @@ public:
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.
* 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.
@ -509,50 +480,45 @@ public:
*
* units = J kmol-1
*
* @param deltaG Output vector of deltaG's for reactions
* Length: m_ii.
* @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 electrochemical free energy change.
//! Return the vector of values for the reaction electrochemical free
//! energy change.
/*!
* These values depend upon the concentration of the solution and
* the voltage of the phases
* These values depend upon the concentration of the solution and the
* voltage of the phases
*
* units = J kmol-1
*
* @param deltaM Output vector of deltaM's for reactions
* Length: m_ii.
* @param deltaM Output vector of deltaM's for reactions Length: m_ii.
*/
virtual void getDeltaElectrochemPotentials(doublereal* deltaM) {
err("getDeltaElectrochemPotentials");
}
/**
* Return the vector of values for the reactions change in
* enthalpy. These values depend upon the concentration of
* the solution.
* 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.
* @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.
* 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.
* @param deltaS Output vector of deltaS's for reactions Length: m_ii.
*/
virtual void getDeltaEntropy(doublereal* deltaS) {
err("getDeltaEntropy");
@ -565,8 +531,7 @@ public:
*
* units = J kmol-1
*
* @param deltaG Output vector of ss deltaG's for reactions
* Length: m_ii.
* @param deltaG Output vector of ss deltaG's for reactions Length: m_ii.
*/
virtual void getDeltaSSGibbs(doublereal* deltaG) {
err("getDeltaSSGibbs");
@ -579,8 +544,7 @@ public:
*
* units = J kmol-1
*
* @param deltaH Output vector of ss deltaH's for reactions
* Length: m_ii.
* @param deltaH Output vector of ss deltaH's for reactions Length: m_ii.
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH) {
err("getDeltaSSEnthalpy");
@ -593,41 +557,33 @@ public:
*
* units = J kmol-1 Kelvin-1
*
* @param deltaS Output vector of ss deltaS's for reactions
* Length: m_ii.
* @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.
* 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.
* @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.
* 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.
* @param ddot Output vector of destruction rates. Length: m_kk.
*/
virtual void getDestructionRates(doublereal* ddot) {
err("getDestructionRates");
@ -639,24 +595,18 @@ public:
* 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.
* @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.
* Stoichiometric coefficient of species k as a reactant in reaction i.
*
* @param k kinetic species index
* @param i reaction index
@ -667,8 +617,7 @@ public:
}
/**
* Stoichiometric coefficient of species k as a product in
* reaction i.
* Stoichiometric coefficient of species k as a product in reaction i.
*
* @param k kinetic species index
* @param i reaction index
@ -708,8 +657,8 @@ public:
//! Get the vector of activity concentrations used in the kinetics object
/*!
* @param conc (output) Vector of activity concentrations. Length is
* equal to the number of species in the kinetics object
* @param[out] conc Vector of activity concentrations. Length is equal
* to the number of species in the kinetics object
*/
virtual void getActivityConcentrations(doublereal* const conc) {
err("getActivityConcentrations");
@ -760,7 +709,7 @@ public:
}
/**
* Return a std::string representing the reaction.
* Return a string representing the reaction.
*
* @param i reaction index
*/
@ -772,12 +721,10 @@ public:
/**
* Return the forward rate constants
*
* length is the number of reactions. units depends
* on many issues. @todo DGG: recommend changing name to
* getFwdRateCoefficients.
* length is the number of reactions. units depends on many issues.
*
* @param kfwd Output vector containing the forward reaction rate constants.
* Length: m_ii.
* @param kfwd Output vector containing the forward reaction rate
* constants. Length: m_ii.
*/
virtual void getFwdRateConstants(doublereal* kfwd) {
err("getFwdRateConstants");
@ -786,11 +733,9 @@ public:
/**
* 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.
* 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.
*
* @param krev Output vector of reverse rate constants.
* @param doIrreversible boolean indicating whether irreversible reactions
@ -801,66 +746,56 @@ public:
err("getFwdRateConstants");
}
/**
* Return the activation energies in Kelvin.
*
* length is the number of reactions
*
* @param E Ouptut vector of activation energies.
* Length: m_ii.
* @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.
* 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.
* Prepare the class for the addition of reactions. This method is called
* by 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.
* Finish adding reactions and prepare for use. This method is called by
* 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();
@ -885,18 +820,15 @@ public:
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.
*/
//@{
@ -936,20 +868,20 @@ public:
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
/*!
* Takes as input an array of properties for all species in the mechanism
* and copies those values belonging 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 selectPhase(const doublereal* data, const thermo_t* phase,
doublereal* phase_data);
protected:
//! Number of reactions in the mechanism
size_t m_ii;
@ -959,7 +891,6 @@ protected:
/// Vector of perturbation factors for each reaction's rate of
/// progress vector. It is initialized to one.
///
vector_fp m_perturb;
/**
@ -988,7 +919,8 @@ protected:
*/
std::vector<std::vector<size_t> > m_products;
//! m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator
//! m_thermo is a vector of pointers to ThermoPhase objects that are
//! involved with this kinetics operator
/*!
* For homogeneous kinetics applications, this vector
* will only have one entry. For interfacial reactions, this
@ -999,9 +931,8 @@ protected:
* 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.
* Note that this kinetics object doesn't own these ThermoPhase objects
* and is not responsible for creating or deleting them.
*/
std::vector<thermo_t*> m_thermo;
@ -1023,9 +954,6 @@ protected:
std::map<std::string, size_t> m_phaseindex;
//! Index in the list of phases of the one surface phase.
/*!
*
*/
size_t m_surfphase;
//! Phase Index where reactions are assumed to be taking place
@ -1039,18 +967,15 @@ protected:
size_t m_mindim;
private:
//! Vector of group lists
std::vector<grouplist_t> m_dummygroups;
//! Private function of the class Kinetics, indicating that a function
//! inherited from the base class hasn't had a definition assigned to it
//! Function indicating that a function inherited from the base class
//! hasn't had a definition assigned to it
/*!
* @param m String message
*/
void err(const std::string& m) const;
};
}

128
src/kinetics/AqueousKinetics.cpp
View file

@ -21,10 +21,7 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
/**
* Construct an empty reaction mechanism.
*/
AqueousKinetics::AqueousKinetics(thermo_t* thermo) :
Kinetics(),
m_nfall(0),
@ -38,7 +35,7 @@ AqueousKinetics::AqueousKinetics(thermo_t* thermo) :
addPhase(*thermo);
}
}
//====================================================================================================================
AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) :
Kinetics(),
m_nfall(0),
@ -50,11 +47,11 @@ AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) :
{
*this = right;
}
//====================================================================================================================
AqueousKinetics::~AqueousKinetics()
{
}
//====================================================================================================================
AqueousKinetics& AqueousKinetics::operator=(const AqueousKinetics& right)
{
if (this == &right) {
@ -100,7 +97,7 @@ AqueousKinetics& AqueousKinetics::operator=(const AqueousKinetics& right)
return *this;
}
//====================================================================================================================
Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
AqueousKinetics* gK = new AqueousKinetics(*this);
@ -108,11 +105,6 @@ Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & t
return gK;
}
//====================================================================================================================
/**
* Update temperature-dependent portions of reaction rates and
* falloff functions.
*/
void AqueousKinetics::
update_T() {}
@ -128,13 +120,8 @@ void AqueousKinetics::_update_rates_T()
m_temp = T;
updateKc();
m_ROP_ok = false;
};
}
/**
* Update properties that depend on concentrations. Currently only
* the enhanced collision partner concentrations are updated here.
*/
void AqueousKinetics::
_update_rates_C()
{
@ -143,9 +130,6 @@ _update_rates_C()
m_ROP_ok = false;
}
/**
* Update the equilibrium constants in molar units.
*/
void AqueousKinetics::updateKc()
{
doublereal rt = GasConstant * m_temp;
@ -172,10 +156,6 @@ void AqueousKinetics::updateKc()
}
}
/**
* Get the equilibrium constants of all reactions, whether
* reversible or not.
*/
void AqueousKinetics::getEquilibriumConstants(doublereal* kc)
{
_update_rates_T();
@ -201,17 +181,6 @@ void AqueousKinetics::getEquilibriumConstants(doublereal* kc)
m_temp = 0.0;
}
/**
*
* getDeltaGibbs():
*
* Return the vector of values for the reaction gibbs free energy
* change
* These values depend upon the concentration
* of the ideal gas.
*
* units = J kmol-1
*/
void AqueousKinetics::getDeltaGibbs(doublereal* deltaG)
{
/*
@ -226,17 +195,6 @@ void AqueousKinetics::getDeltaGibbs(doublereal* deltaG)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
/**
*
* 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
*/
void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
/*
@ -251,17 +209,6 @@ void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
/*
*
* 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
*/
void AqueousKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
@ -276,17 +223,6 @@ void AqueousKinetics::getDeltaEntropy(doublereal* deltaS)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
/**
*
* 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
*/
void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
/*
@ -303,17 +239,6 @@ void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
/**
*
* 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
*/
void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
@ -334,17 +259,6 @@ void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
/*
*
* 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
*/
void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
@ -364,8 +278,6 @@ void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void AqueousKinetics::updateROP()
{
_update_rates_T();
@ -405,14 +317,6 @@ void AqueousKinetics::updateROP()
m_ROP_ok = true;
}
/**
*
* getFwdRateConstants():
*
* Update the rate of progress for the reactions.
* This key routine makes sure that the rate of progress vectors
* located in the solid kinetics data class are up to date.
*/
void AqueousKinetics::
getFwdRateConstants(doublereal* kfwd)
{
@ -430,17 +334,6 @@ getFwdRateConstants(doublereal* kfwd)
}
}
/**
*
* getRevRateConstants():
*
* Return a vector of the reverse reaction 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.
*/
void AqueousKinetics::
getRevRateConstants(doublereal* krev, bool doIrreversible)
{
@ -468,7 +361,6 @@ getRevRateConstants(doublereal* krev, bool doIrreversible)
void AqueousKinetics::addReaction(ReactionData& r)
{
if (r.reactionType == ELEMENTARY_RXN) {
addElementaryReaction(r);
}
@ -480,9 +372,6 @@ void AqueousKinetics::addReaction(ReactionData& r)
m_rxneqn.push_back(r.equation);
}
void AqueousKinetics::addElementaryReaction(ReactionData& r)
{
size_t iloc;
@ -498,9 +387,6 @@ void AqueousKinetics::addElementaryReaction(ReactionData& r)
registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc);
}
void AqueousKinetics::installReagents(const ReactionData& r)
{
@ -568,7 +454,6 @@ void AqueousKinetics::installReagents(const ReactionData& r)
}
}
void AqueousKinetics::installGroups(size_t irxn,
const vector<grouplist_t>& r,
const vector<grouplist_t>& p)
@ -580,7 +465,6 @@ void AqueousKinetics::installGroups(size_t irxn,
}
}
void AqueousKinetics::init()
{
m_kk = thermo().nSpecies();

23
src/kinetics/GRI_30_Kinetics.cpp
View file

@ -21,9 +21,6 @@ using namespace std;
namespace Cantera
{
/**
* Construct an empty reaction mechanism.
*/
GRI_30_Kinetics::
GRI_30_Kinetics(thermo_t* th) : GasKinetics(th) {}
@ -40,13 +37,8 @@ gri30_update_rates_T()
m_temp = T;
gri30_updateKc();
m_ROP_ok = false;
//}
};
}
/**
* Update the equilibrium constants in molar units.
*/
void GRI_30_Kinetics::gri30_updateKc()
{
vector_fp a(m_kk);
@ -58,7 +50,6 @@ void GRI_30_Kinetics::gri30_updateKc()
update_kc(&a[0], exp_c_ref, &m_rkcn[0]);
}
void GRI_30_Kinetics::gri30_updateROP()
{
@ -82,7 +73,6 @@ void GRI_30_Kinetics::gri30_updateROP()
m_ROP_ok = true;
}
void GRI_30_Kinetics::update_rates(doublereal t, doublereal tlog, doublereal* rf)
{
doublereal rt = 1.0/t;
@ -289,7 +279,6 @@ void GRI_30_Kinetics::update_rates(doublereal t, doublereal tlog, doublereal* rf
rf[324] = exp(23.6818 + -0.32 * tlog);
}
void GRI_30_Kinetics::update_kc(const doublereal* a, doublereal exp_c0, doublereal* rkc)
{
rkc[0] = a[3]*exp_c0/(a[2]*a[2]);
@ -603,7 +592,6 @@ void GRI_30_Kinetics::update_kc(const doublereal* a, doublereal exp_c0, doubler
rkc[324] = a[25]*a[25]/(a[12]*a[49]);
}
void GRI_30_Kinetics::get_wdot(const doublereal* rop, doublereal* wdot)
{
wdot[0] = - rop[2] + rop[7] + rop[38] + rop[39] + rop[40] + rop[41] + rop[44] + rop[46] + rop[48] + rop[50] + rop[52] + rop[54] + rop[57] + rop[59] + rop[64] + rop[67] + rop[68] + rop[72] + rop[74] + rop[76] + rop[77] + rop[79] - rop[82] - rop[83] - rop[125] - rop[135] + rop[136] - rop[145] - rop[171] + rop[173] + rop[190] + rop[196] + rop[201] + rop[208] + rop[213] - rop[220] + rop[265] + rop[275] + rop[276] + rop[283] + rop[287] - rop[288] + rop[292] + rop[298] + rop[299] + rop[308] + rop[313];
@ -661,7 +649,6 @@ void GRI_30_Kinetics::get_wdot(const doublereal* rop, doublereal* wdot)
wdot[52] = + rop[285] - rop[295] - rop[296] - rop[297] - rop[298] - rop[299] - rop[300] - rop[301] - rop[302];
}
void GRI_30_Kinetics::eval_ropnet(const doublereal* c, const doublereal* rf, const doublereal* rkc, doublereal* r)
{
r[0] = rf[0] * (c[2] * c[2] - rkc[0] * c[3]);
@ -992,11 +979,3 @@ void GRI_30_Kinetics::eval_ropnet(const doublereal* c, const doublereal* rf, con
}
}

199
src/kinetics/GasKinetics.cpp
View file

@ -2,7 +2,6 @@
* @file GasKinetics.cpp
*
* Homogeneous kinetics in ideal gases
*
*/
// Copyright 2001 California Institute of Technology
@ -19,11 +18,6 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
/*
* Construct an empty reaction mechanism.
*/
GasKinetics::
GasKinetics(thermo_t* thermo) :
Kinetics(),
@ -43,7 +37,6 @@ GasKinetics(thermo_t* thermo) :
m_temp = 0.0;
}
//====================================================================================================================
GasKinetics::GasKinetics(const GasKinetics& right) :
Kinetics(),
m_nfall(0),
@ -59,11 +52,11 @@ GasKinetics::GasKinetics(const GasKinetics& right) :
m_temp = 0.0;
*this = right;
}
//====================================================================================================================
GasKinetics::~GasKinetics()
{
}
//====================================================================================================================
GasKinetics& GasKinetics::operator=(const GasKinetics& right)
{
if (this == &right) {
@ -124,26 +117,14 @@ GasKinetics& GasKinetics::operator=(const GasKinetics& right)
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.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
*/
Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
GasKinetics* gK = new GasKinetics(*this);
gK->assignShallowPointers(tpVector);
return gK;
}
//====================================================================================================================
void GasKinetics::update_rates_T()
{
doublereal T = thermo().temperature();
@ -171,9 +152,7 @@ void GasKinetics::update_rates_T()
m_temp = T;
updateKc();
m_ROP_ok = false;
};
//====================================================================================================================
}
void GasKinetics::update_rates_C()
{
@ -204,10 +183,7 @@ void GasKinetics::update_rates_C()
m_ROP_ok = false;
}
//====================================================================================================================
/**
* Update the equilibrium constants in molar units.
*/
void GasKinetics::updateKc()
{
thermo().getStandardChemPotentials(&m_grt[0]);
@ -227,11 +203,7 @@ void GasKinetics::updateKc()
m_rkcn[ m_irrev[i] ] = 0.0;
}
}
//====================================================================================================================
/**
* Get the equilibrium constants of all reactions, whether
* reversible or not.
*/
void GasKinetics::getEquilibriumConstants(doublereal* kc)
{
update_rates_T();
@ -250,18 +222,7 @@ void GasKinetics::getEquilibriumConstants(doublereal* kc)
// be updated before it is used next.
m_temp = 0.0;
}
//====================================================================================================================
/**
*
* getDeltaGibbs():
*
* Return the vector of values for the reaction gibbs free energy
* change
* These values depend upon the concentration
* of the ideal gas.
*
* units = J kmol-1
*/
void GasKinetics::getDeltaGibbs(doublereal* deltaG)
{
/*
@ -275,18 +236,7 @@ void GasKinetics::getDeltaGibbs(doublereal* deltaG)
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
//====================================================================================================================
/**
*
* 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
*/
void GasKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
/*
@ -300,18 +250,7 @@ void GasKinetics::getDeltaEnthalpy(doublereal* deltaH)
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
//====================================================================================================================
/*
*
* 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
*/
void GasKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
@ -325,18 +264,7 @@ void GasKinetics::getDeltaEntropy(doublereal* deltaS)
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
//====================================================================================================================
/**
*
* 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
*/
void GasKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
/*
@ -352,18 +280,7 @@ void GasKinetics::getDeltaSSGibbs(doublereal* deltaG)
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
//====================================================================================================================
/**
*
* 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
*/
void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
@ -383,18 +300,7 @@ void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
//====================================================================================================================
/*********************************************************************
*
* 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
*/
void GasKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
@ -414,57 +320,24 @@ void GasKinetics::getDeltaSSEntropy(doublereal* deltaS)
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
//====================================================================================================================
// Return the species net production rates
/*
* Species net production rates [kmol/m^3/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.
*
* @param net Array of species production rates.
* units kmol m-3 s-1
*/
void GasKinetics::getNetProductionRates(doublereal* net)
{
updateROP();
m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net);
}
//====================================================================================================================
// Return the species creation rates
/*
* Species creation rates [kmol/m^3]. Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species.
*
* @param cdot Array of species production rates.
* units kmol m-3 s-1
*/
void GasKinetics::getCreationRates(doublereal* cdot)
{
updateROP();
m_rxnstoich.getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot);
}
//====================================================================================================================
// Return a vector of the species destruction rates
/*
* Species destruction rates [kmol/m^3]. Return the species
* destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species.
*
*
* @param ddot Array of species destruction rates.
* units kmol m-3 s-1
*
*/
void GasKinetics::getDestructionRates(doublereal* ddot)
{
updateROP();
m_rxnstoich.getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot);
}
//====================================================================================================================
void GasKinetics::processFalloffReactions()
{
// use m_ropr for temporary storage of reduced pressure
@ -485,7 +358,6 @@ void GasKinetics::processFalloffReactions()
m_ropf.begin(), m_fallindx.begin());
}
//====================================================================================================================
void GasKinetics::updateROP()
{
update_rates_C();
@ -533,15 +405,7 @@ void GasKinetics::updateROP()
m_ROP_ok = true;
}
//====================================================================================================================
/**
*
* getFwdRateConstants():
*
* Update the rate of progress for the reactions.
* This key routine makes sure that the rate of progress vectors
* located in the solid kinetics data class are up to date.
*/
void GasKinetics::
getFwdRateConstants(doublereal* kfwd)
{
@ -571,18 +435,7 @@ getFwdRateConstants(doublereal* kfwd)
kfwd[i] = m_ropf[i];
}
}
//====================================================================================================================
/**
*
* getRevRateConstants():
*
* Return a vector of the reverse reaction 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.
*/
void GasKinetics::
getRevRateConstants(doublereal* krev, bool doIrreversible)
{
@ -605,7 +458,7 @@ getRevRateConstants(doublereal* krev, bool doIrreversible)
}
}
}
//====================================================================================================================
void GasKinetics::
addReaction(ReactionData& r)
{
@ -636,7 +489,6 @@ addReaction(ReactionData& r)
m_rxneqn.push_back(r.equation);
}
//====================================================================================================================
void GasKinetics::
addFalloffReaction(ReactionData& r)
{
@ -673,7 +525,6 @@ addFalloffReaction(ReactionData& r)
++m_nfall;
registerReaction(reactionNumber(), FALLOFF_RXN, iloc);
}
//====================================================================================================================
void GasKinetics::
addElementaryReaction(ReactionData& r)
@ -689,7 +540,6 @@ addElementaryReaction(ReactionData& r)
registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc);
}
//====================================================================================================================
void GasKinetics::
addThreeBodyReaction(ReactionData& r)
{
@ -706,7 +556,6 @@ addThreeBodyReaction(ReactionData& r)
r.default_3b_eff);
registerReaction(reactionNumber(), THREE_BODY_RXN, iloc);
}
//====================================================================================================================
void GasKinetics::addPlogReaction(ReactionData& r)
{
@ -795,8 +644,6 @@ void GasKinetics::installReagents(const ReactionData& r)
m_nirrev++;
}
}
//====================================================================================================================
void GasKinetics::installGroups(size_t irxn,
const vector<grouplist_t>& r, const vector<grouplist_t>& p)
@ -808,7 +655,6 @@ void GasKinetics::installGroups(size_t irxn,
}
}
//====================================================================================================================
void GasKinetics::init()
{
m_kk = thermo().nSpecies();
@ -818,7 +664,7 @@ void GasKinetics::init()
m_grt.resize(m_kk);
m_logp_ref = log(thermo().refPressure()) - log(GasConstant);
}
//====================================================================================================================
void GasKinetics::finalize()
{
if (!m_finalized) {
@ -838,11 +684,10 @@ void GasKinetics::finalize()
}
}
}
//====================================================================================================================
bool GasKinetics::ready() const
{
return m_finalized;
}
//====================================================================================================================
}
//======================================================================================================================

328
src/kinetics/InterfaceKinetics.cpp
View file

@ -1,6 +1,5 @@
/**
* @file InterfaceKinetics.cpp
*
*/
// Copyright 2002 California Institute of Technology
@ -18,16 +17,6 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
/*
* Construct an empty InterfaceKinetics reaction mechanism.
* @param thermo This is an optional parameter that may be
* used to initialize the inherited Kinetics class with
* one ThermoPhase class object -> in other words it's
* useful for initialization of homogeneous kinetics
* mechanisms.
*/
InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
Kinetics(),
m_redo_rates(false),
@ -61,20 +50,12 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
addPhase(*thermo);
}
}
//====================================================================================================================
/*
* Destructor
*/
InterfaceKinetics::~InterfaceKinetics()
{
delete m_integrator;
}
//====================================================================================================================
// 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_redo_rates(false),
@ -109,14 +90,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
*/
*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)
{
@ -183,48 +157,20 @@ int InterfaceKinetics::type() const
{
return cInterfaceKinetics;
}
//====================================================================================================================
// 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.
*
* @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
* m_thermo vector within this object
*/
Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
InterfaceKinetics* iK = new InterfaceKinetics(*this);
iK->assignShallowPointers(tpVector);
return iK;
}
//====================================================================================================================
// Set the electric potential in the nth phase
/*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void InterfaceKinetics::setElectricPotential(int n, doublereal V)
{
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
//====================================================================================================================
// Update properties that depend on temperature
/*
* This is called to update all of the properties that depend on temperature
*
* Current objects that this function updates
* m_logtemp
* m_rfn
* m_rates.
* updateKc();
*/
void InterfaceKinetics::_update_rates_T()
{
_update_rates_phi();
@ -250,7 +196,7 @@ void InterfaceKinetics::_update_rates_T()
m_redo_rates = false;
}
}
//====================================================================================================================
void InterfaceKinetics::_update_rates_phi()
{
for (size_t n = 0; n < nPhases(); n++) {
@ -260,16 +206,7 @@ void InterfaceKinetics::_update_rates_phi()
}
}
}
//====================================================================================================================
/**
* Update properties that depend on concentrations. This method
* fills out the array of generalized concentrations by calling
* method getActivityConcentrations for each phase, which classes
* representing phases should overload to return the appropriate
* quantities.
*/
void InterfaceKinetics::_update_rates_C()
{
for (size_t n = 0; n < nPhases(); n++) {
@ -286,26 +223,12 @@ void InterfaceKinetics::_update_rates_C()
m_ROP_ok = false;
}
// Get the vector of activity concentrations used in the kinetics object
/*
* @param conc (output) Vector of activity concentrations. Length is
* equal to the number of species in the kinetics object
*/
void InterfaceKinetics::getActivityConcentrations(doublereal* const conc)
{
_update_rates_C();
copy(m_conc.begin(), m_conc.end(), conc);
}
/**
* Update the equilibrium constants in molar units for all
* reversible reactions. Irreversible reactions have their
* equilibrium constant set to zero.
* For reactions involving charged species the equilibrium
* constant is adjusted according to the electrostatic potential.
*/
void InterfaceKinetics::updateKc()
{
fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
@ -347,9 +270,6 @@ void InterfaceKinetics::updateKc()
}
}
}
//====================================================================================================================
void InterfaceKinetics::checkPartialEquil()
{
@ -405,10 +325,6 @@ void InterfaceKinetics::getNetRatesOfProgress(doublereal* netROP)
std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
}
/**
* Get the equilibrium constants of all reactions, whether
* reversible or not.
*/
void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
{
size_t ik=0;
@ -462,62 +378,24 @@ void InterfaceKinetics::getExchangeCurrentQuantities()
}
// Returns the Species creation rates [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 of the kinetics
* model
*
* @param cdot Vector containing the creation rates.
* length = m_kk. units = kmol/m^2/s
*/
void InterfaceKinetics::getCreationRates(doublereal* cdot)
{
updateROP();
m_rxnstoich.getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot);
}
// Return the Species destruction rates [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 in all phases of the kinetics model
*/
void InterfaceKinetics::getDestructionRates(doublereal* ddot)
{
updateROP();
m_rxnstoich.getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot);
}
// Return the species net production rates
/*
* Species net production rates [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 in all phases of the kinetics
* model
*
* @param net Vector of species production rates.
* units kmol m-d s-1, where d is dimension.
*/
void InterfaceKinetics::getNetProductionRates(doublereal* net)
{
updateROP();
m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net);
}
//====================================================================================================================
// Apply corrections for interfacial charge transfer reactions
/*
* For reactions that transfer charge across a potential difference,
* the activation energies are modified by the potential difference.
* (see, for example, ...). This method applies this correction.
*
* @param kf Vector of forward reaction rate constants on which to have
* the correction applied
*/
void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf)
{
// compute the electrical potential energy of each species
@ -573,7 +451,7 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf)
}
}
}
//====================================================================================================================
void InterfaceKinetics::applyExchangeCurrentDensityFormulation(doublereal* const kfwd)
{
getExchangeCurrentQuantities();
@ -589,14 +467,8 @@ void InterfaceKinetics::applyExchangeCurrentDensityFormulation(doublereal* const
kfwd[irxn] *= tmp;
}
}
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaction
* mechanism. This routine operates on internal data.
*/
void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd)
{
@ -609,12 +481,7 @@ void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd)
multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin());
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaction
* mechanism. This routine operates on internal data.
*/
void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
{
getFwdRateConstants(krev);
@ -627,17 +494,12 @@ void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversibl
multiply_each(krev, krev + nReactions(), m_rkcn.begin());
}
}
//====================================================================================================================
void InterfaceKinetics::getActivationEnergies(doublereal* E)
{
copy(m_E.begin(), m_E.end(), E);
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaction
* mechanism. This routine operates on internal data.
*/
void InterfaceKinetics::updateROP()
{
_update_rates_T();
@ -768,19 +630,7 @@ InterfaceKinetics::adjustRatesForIntermediatePhases()
}
#endif
//=================================================================================================
/*
*
* getDeltaGibbs():
*
* Return the vector of values for the reaction gibbs free energy
* change
* These values depend upon the concentration
* of the ideal gas.
*
* units = J kmol-1
*/
void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
{
/*
@ -800,17 +650,7 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
}
//=================================================================================================
// Return the vector of values for the reaction electrochemical free energy change.
/*
* These values depend upon the concentration of the solution and
* the voltage of the phases
*
* units = J kmol-1
*
* @param deltaM Output vector of deltaM's for reactions
* Length: m_ii.
*/
void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
{
/*
@ -827,18 +667,7 @@ void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaM);
}
//=================================================================================================
/*
*
* 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
*/
void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
/*
@ -855,19 +684,6 @@ void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
}
// Return the vector of values for the change in
// entropy due to each reaction
/*
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*
* @param deltaS vector of Enthalpy changes
* Length = m_ii, number of reactions
*
*/
void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
@ -884,17 +700,6 @@ void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
}
/**
*
* 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
*/
void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
/*
@ -913,17 +718,6 @@ void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaG)
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
}
/**
*
* 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
*/
void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
@ -946,17 +740,6 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
}
/*********************************************************************
*
* 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
*/
void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
@ -978,23 +761,8 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
}
//====================================================================================================================
/**
* Add a single reaction to the mechanism. This routine
* must be called after init() and before finalize().
* This function branches on the types of reactions allowed
* by the interfaceKinetics manager in order to install
* the reaction correctly in the manager.
* The manager allows the following reaction types
* Elementary
* Surface
* Global
* There is no difference between elementary and surface
* reactions.
*/
void InterfaceKinetics::addReaction(ReactionData& r)
{
/*
* Install the rate coefficient for the current reaction
* in the appropriate data structure.
@ -1036,7 +804,7 @@ void InterfaceKinetics::addReaction(ReactionData& r)
m_rxnPhaseIsProduct[i][p] = true;
}
}
//====================================================================================================================
void InterfaceKinetics::addElementaryReaction(ReactionData& r)
{
// install rate coeff calculator
@ -1089,7 +857,6 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& r)
m_rfn.push_back(r.rateCoeffParameters[0]);
registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc);
}
//====================================================================================================================
void InterfaceKinetics::setIOFlag(int ioFlag)
{
@ -1126,7 +893,6 @@ void InterfaceKinetics::setIOFlag(int ioFlag)
// registerReaction( reactionNumber(), GLOBAL_RXN, iloc);
// }
void InterfaceKinetics::installReagents(const ReactionData& r)
{
@ -1224,22 +990,14 @@ void InterfaceKinetics::installReagents(const ReactionData& r)
m_nirrev++;
}
}
//===============================================================================================
void InterfaceKinetics::addPhase(thermo_t& thermo)
{
Kinetics::addPhase(thermo);
m_phaseExists.push_back(true);
m_phaseIsStable.push_back(true);
}
//================================================================================================
/**
* Prepare the class for the addition of reactions. This function
* must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
void InterfaceKinetics::init()
{
m_kk = 0;
@ -1254,15 +1012,7 @@ void InterfaceKinetics::init()
m_pot.resize(m_kk, 0.0);
m_phi.resize(nPhases(), 0.0);
}
//================================================================================================
/**
* Finish adding reactions and prepare for use. This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*
* Here, we resize work arrays based on the number of reactions,
* since we don't know this number up to now.
*/
void InterfaceKinetics::finalize()
{
Kinetics::finalize();
@ -1308,16 +1058,11 @@ doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
return 0.0;
}
//================================================================================================
bool InterfaceKinetics::ready() const
{
return m_finalized;
}
//================================================================================================
// Advance the surface coverages in time
/*
* @param tstep Time value to advance the surface coverages
*/
void InterfaceKinetics::
advanceCoverages(doublereal tstep)
{
@ -1331,19 +1076,7 @@ advanceCoverages(doublereal tstep)
delete m_integrator;
m_integrator = 0;
}
//================================================================================================
// Solve for the pseudo steady-state of the surface problem
/*
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*
* the integrator object is saved between calls to
* reduce the computational cost of repeated calls.
*/
void InterfaceKinetics::
solvePseudoSteadyStateProblem(int ifuncOverride, doublereal timeScaleOverride)
{
@ -1360,7 +1093,6 @@ solvePseudoSteadyStateProblem(int ifuncOverride, doublereal timeScaleOverride)
*/
m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
}
//================================================================================================
void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
{
@ -1385,14 +1117,7 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
}
}
//================================================================================================
// Gets the phase existence int for the ith phase
/*
* @param iphase Phase Id
*
* @return Returns the int specifying whether the kinetics object thinks the phase exists
* or not. If it exists, then species in that phase can be a reactant in reactions.
*/
int InterfaceKinetics::phaseExistence(const size_t iphase) const
{
if (iphase >= m_thermo.size()) {
@ -1400,16 +1125,7 @@ int InterfaceKinetics::phaseExistence(const size_t iphase) const
}
return m_phaseExists[iphase];
}
//================================================================================================
// Gets the phase stability int for the ith phase
/*
* @param iphase Phase Id
*
* @return Returns the int specifying whether the kinetics object thinks the phase is stable
* with nonzero mole numbers.
* If it stable, then the kinetics object will allow for rates of production of
* of species in that phase that are positive.
*/
int InterfaceKinetics::phaseStability(const size_t iphase) const
{
if (iphase >= m_thermo.size()) {
@ -1417,7 +1133,6 @@ int InterfaceKinetics::phaseStability(const size_t iphase) const
}
return m_phaseIsStable[iphase];
}
//================================================================================================
void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
{
@ -1431,7 +1146,6 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
}
}
//================================================================================================
void EdgeKinetics::finalize()
{
m_rwork.resize(std::max<size_t>(nReactions(), 1));
@ -1456,7 +1170,5 @@ void EdgeKinetics::finalize()
m_finalized = true;
}
//================================================================================================
}

112
src/kinetics/Kinetics.cpp
View file

@ -8,8 +8,6 @@
*/
// Copyright 2001-2004 California Institute of Technology
// Why InterfaceKinetics.h and not Kinetics.h ??
#include "cantera/kinetics/Kinetics.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/kinetics/StoichManager.h"
@ -23,8 +21,6 @@ using namespace std;
namespace Cantera
{
Kinetics::Kinetics() :
m_ii(0),
m_kk(0),
@ -43,12 +39,6 @@ Kinetics::Kinetics() :
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_kk(0),
@ -69,13 +59,6 @@ Kinetics::Kinetics(const Kinetics& right) :
*this = 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)
{
@ -104,17 +87,6 @@ operator=(const Kinetics& right)
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 std::vector<thermo_t*> & tpVector) const
{
Kinetics* ko = new Kinetics(*this);
@ -170,7 +142,6 @@ void Kinetics::checkSpeciesArraySize(size_t kk) const
}
}
//====================================================================================================================
void Kinetics::assignShallowPointers(const std::vector<thermo_t*> & tpVector)
{
size_t ns = tpVector.size();
@ -198,17 +169,7 @@ void Kinetics::assignShallowPointers(const std::vector<thermo_t*> & tpVector)
}
//====================================================================================================================
/**
* Takes as input an array of properties for all species in the
* mechanism and copies those values belonging 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)
{
@ -223,17 +184,6 @@ void Kinetics::selectPhase(const doublereal* data, const thermo_t* phase,
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
* "<unknown>" is returned.
*/
string Kinetics::kineticsSpeciesName(size_t k) const
{
for (size_t n = m_start.size()-1; n != npos; n--) {
@ -244,17 +194,6 @@ string Kinetics::kineticsSpeciesName(size_t k) const
return "<unknown>";
}
/**
* 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.
*
* return
* - If a match is found, the position in the species list is returned.
* - If no match is found, the value -1 is returned.
*
* @param nm Input string name of the species
*/
size_t Kinetics::kineticsSpeciesIndex(const std::string& nm) const
{
for (size_t n = 0; n < m_thermo.size(); n++) {
@ -268,18 +207,6 @@ size_t Kinetics::kineticsSpeciesIndex(const std::string& nm) const
return npos;
}
/**
* This routine will look up a species number based on the input
* std::string nm. The lookup of species will occur in the specified
* phase of the object, or all phases if ph is "<any>".
*
* return
* - If a match is found, the position in the species list is returned.
* - If no match is found, the value npos (-1) is returned.
*
* @param nm Input string name of the species
* @param ph Input string name of the phase.
*/
size_t Kinetics::kineticsSpeciesIndex(const std::string& nm,
const std::string& ph) const
{
@ -300,13 +227,6 @@ size_t Kinetics::kineticsSpeciesIndex(const std::string& nm,
return npos;
}
/**
* 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(const std::string& nm)
{
size_t np = m_thermo.size();
@ -322,13 +242,6 @@ thermo_t& Kinetics::speciesPhase(const std::string& nm)
return thermo(0);
}
//==============================================================================================
/*
* 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.
*/
size_t Kinetics::speciesPhaseIndex(size_t k)
{
for (size_t n = m_start.size()-1; n != npos; n--) {
@ -340,26 +253,8 @@ size_t Kinetics::speciesPhaseIndex(size_t k)
return npos;
}
/*
* 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) {
@ -402,11 +297,6 @@ void Kinetics::finalize()
}
}
// Private function of the class Kinetics, indicating that a function
// inherited from the base class hasn't had a definition assigned to it
/*
* @param m String message
*/
void Kinetics::err(const std::string& m) const
{
throw CanteraError("Kinetics::" + m,