diff --git a/include/cantera/kinetics/ElectrodeKinetics.h b/include/cantera/kinetics/ElectrodeKinetics.h new file mode 100644 index 000000000..f49530d0d --- /dev/null +++ b/include/cantera/kinetics/ElectrodeKinetics.h @@ -0,0 +1,87 @@ +/** + * @file ElectrodeKinetics.h + * + * @ingroup chemkinetics + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef CT_ELECTRODEKINETICS_H +#define CT_ELECTRODEKINETICS_H + +#include "InterfaceKinetics.h" + + +namespace Cantera +{ + + +//! A kinetics manager for heterogeneous reaction mechanisms. The +//! reactions are assumed to occur at a 2D interface between two 3D phases. +/*! + * This class is a slight addition to the InterfaceKinetics class, adding + * several concepts. First we explicity identify the electrode and solution + * phases. We will also assume that there is an electron phase. + * + * @ingroup chemkinetics + */ +class ElectrodeKinetics : public InterfaceKinetics +{ +public: + //! Constructor + /*! + * @param thermo The optional parameter may be used to initialize + * the object with one ThermoPhase object. + * HKM Note -> Since the interface kinetics + * object will probably require multiple thermophase + * objects, this is probably not a good idea + * to have this parameter. + */ + ElectrodeKinetics(thermo_t* thermo = 0); + + /// Destructor. + virtual ~ElectrodeKinetics(); + + //! Copy Constructor for the %Kinetics object. + ElectrodeKinetics(const ElectrodeKinetics& right); + + //! Assignment operator + ElectrodeKinetics& operator=(const ElectrodeKinetics& right); + + virtual Kinetics* duplMyselfAsKinetics(const std::vector & tpVector) const; + + virtual int type() const; + + + + + + //void addGlobalReaction(ReactionData& r); + + + + + +protected: + + //! index of the metal phase in the list of phases for this surface + size_t metalPhaseRS_; + + size_t electronPhaseRS_; + + //! Index of the solution phase in the list of phases for this surface + size_t solnPhaseRS_; + + //! Index of the electrons species in the list of species for this surface kinetics, if none set it to -1 + size_t kElectronRS_; + + + + +}; +} + +#endif diff --git a/include/cantera/kinetics/ExtraGlobalRxn.h b/include/cantera/kinetics/ExtraGlobalRxn.h new file mode 100644 index 000000000..d34feb263 --- /dev/null +++ b/include/cantera/kinetics/ExtraGlobalRxn.h @@ -0,0 +1,120 @@ +/** + * @file ReactingVolDomain.h + * + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef EXTRAGLOBALRXN_H +#define EXTRAGLOBALRXN_H + +#include "cantera/kinetics/InterfaceKinetics.h" +#include +#include + +namespace Cantera +{ + +//! Class describing an extra global reaction, which is defined as +//! a linear combination of actuals reactions, global or mass-action, creating a global stoichiometric result +/*! + * This is useful for defining thermodynamics of global processes that occur + * on a surface or in a homogeneous phase. + * + * The class is set up via the function setupElemRxnVector(RxnVector, specialSpecies) which defines + * the vector of stoichiometric coefficients representing the base reaction to combine in order to + * achieve the global result that's to be calculated. specialSpecies is the index of the species + * within the kinetics object that is used to identify the global reaction. Rates of progress + * are defined in terms of the production rate of the special species. + * + */ +class ExtraGlobalRxn +{ + +public: + //! Constructor takes a default kinetics pointer + /*! + * @param[in] k_ptr Pointer to a Kinetics class that will be used as the basis + * for constructing this class. + */ + ExtraGlobalRxn(Kinetics* k_ptr); + + //! Destructor + virtual ~ExtraGlobalRxn(); + + void setupElemRxnVector(double* RxnVector, + int specialSpecies = -1); + std::string reactionString(); + double deltaSpecValue(double* speciesVectorProperty); + + std::vector& reactants(); + std::vector& products(); + bool isReversible(); + + double ROPValue(double* ROPKinVector); + double FwdROPValue(double* FwdROPElemKinVector, double* RevROPElemKinVector); + double RevROPValue(double* FwdROPElemKinVector, double* RevROPElemKinVector); + + double reactantStoichCoeff(int kKin); + double productStoichCoeff(int kKin); + bool m_ThisIsASurfaceRxn; + double deltaRxnVecValue(double* rxnVectorProperty); + + //! This kinetics operator is associated with just one + //! homogeneous phase, associated with tpList[0] phase + /*! + * Kinetics object pointer + */ + Cantera::Kinetics* m_kinetics; + + //! This kinetics operator is associated with multiple + //! homogeneous and surface phases. + /*! + * This object owns the Kinetics object + */ + Cantera::InterfaceKinetics* m_InterfaceKinetics; + + int m_nKinSpecies; + + //! Number of reactants in the global reaction + int m_nReactants; + + //! Vector of reactants that make up the global reaction + /*! + * This is a list of reactants using the kinetic species index + */ + std::vector m_Reactants; + + //! Vector of reactant stoichiometries that make up the global reaction + /*! + * This is a list of reactant stoichiometries. The species index is given in + * the member m_Reactants using the kinetic species index. + */ + std::vector m_ReactantStoich; + + int m_nProducts; + std::vector m_Products; + std::vector m_ProductStoich; + + int m_nNetSpecies; + std::vector m_NetSpecies; + std::vector m_netStoich; + + int m_nRxns; + std::vector m_ElemRxnVector; + + int m_SpecialSpecies; + bool m_SpecialSpeciesProduct; + int m_SS_index; + + int iphaseKin; + bool m_ok; + bool m_reversible; + + +}; +} +#endif diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index ba6f54dd8..3a7c66f5b 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -105,6 +105,16 @@ public: */ void updateExchangeCurrentQuantities(); + //! Return the vector of values for the reaction gibbs free energy change. + /*! + * (virtual from Kinetics.h) + * 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. + * If 0, this updates the internally stored values only. + */ virtual void getDeltaGibbs(doublereal* deltaG); virtual void getDeltaElectrochemPotentials(doublereal* deltaM); @@ -306,11 +316,13 @@ public: /*! * @param rxnNumber reaction number * @param type reaction type - * @param loc location ?? + * @param loc location location in the reaction rate coefficient calculator + * + * Right now we only use one reaction rate coefficient calculated named ELEMENTARY_RXN + * Therefore, this is not used within the code) + * (type, loc) is stored as a std::pair */ - void registerReaction(size_t rxnNumber, int type, size_t loc) { - m_index[rxnNumber] = std::pair(type, loc); - } + void registerReaction(size_t rxnNumber, int type, size_t loc); //! Apply modifications for the fowward reaction rate for interfacial charge transfer reactions /*! @@ -396,6 +408,8 @@ 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 @@ -515,6 +529,14 @@ protected: */ vector_fp m_mu0; + //! Vector of chemical potentials for all species + /*! + * This vector contains a vector of chemical potentials for all of the species in the kinetics object + * + * Length = m_kk. Units = J/kmol. + */ + vector_fp m_mu; + //! Vector of standard state electrochemical potentials modified by //! a standard concentration term. /*! @@ -636,6 +658,19 @@ protected: * units are Joule kmol-1 */ vector_fp m_deltaG0; + + //! Vector of deltaG[] of reaction, the delta gibbs free energies for each reaction + /*! + * Length is the number of reactions + * units are Joule kmol-1 + */ + vector_fp m_deltaG; + + //! Vector of the products of the standard concentrations of the reactants + /*! + * Units vary wrt what the units of the standard concentrations are + * Length = number of reactions. + */ vector_fp m_ProdStanConcReac; doublereal m_logp0; diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index af7941767..f820331c6 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -465,6 +465,7 @@ public: //! Return the vector of values for the reaction gibbs free energy change. /*! + * (virtual from Kinetics.h) * These values depend upon the concentration of the solution. * * units = J kmol-1 diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index 28587c2d3..dd767cca1 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -13,6 +13,9 @@ namespace Cantera //! Intermediate class which stores data about a reaction and its rate //! parameterization before adding the reaction to a Kinetics object. +/*! + * All data in this class is public. + */ class ReactionData { public: diff --git a/include/cantera/kinetics/ReactionStoichMgr.h b/include/cantera/kinetics/ReactionStoichMgr.h index efd66831b..88168c2b0 100644 --- a/include/cantera/kinetics/ReactionStoichMgr.h +++ b/include/cantera/kinetics/ReactionStoichMgr.h @@ -66,8 +66,9 @@ public: ReactionStoichMgr& operator=(const ReactionStoichMgr& right); - /** - * Add a reaction with mass-action kinetics. Vectors + + //! Add a reaction with mass-action kinetics. + /*!Vectors * 'reactants' and 'products' contain the integer species * indices of the reactants and products, respectively. Note * that if more than one molecule of a given species is @@ -86,11 +87,11 @@ public: * - reactants: (2, 2) [ note repeated index ] * - products: (1) * - * @param rxn Reaction number. This number will be used as the index - * into the rate of progress vector in the methods below. - * @param reactants vector of integer reactant indices - * @param products vector of integer product indices - * @param reversible true if the reaction is reversible, false otherwise + * @param rxn Reaction number. This number will be used as the index + * into the rate of progress vector in the methods below. + * @param reactants Vector of integer reactant indices + * @param products Vector of integer product indices + * @param reversible True if the reaction is reversible, false otherwise */ virtual void add(size_t rxn, const std::vector& reactants, const std::vector& products, bool reversible); diff --git a/include/cantera/kinetics/RxnMolChange.h b/include/cantera/kinetics/RxnMolChange.h new file mode 100644 index 000000000..f842aba4c --- /dev/null +++ b/include/cantera/kinetics/RxnMolChange.h @@ -0,0 +1,117 @@ +/** + * @file RxnMolChange.cpp + * + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef RXNMOLCHANGE_H +#define RXNMOLCHANGE_H + + +#include + +namespace Cantera +{ +class ExtraGlobalRxn; +class Kinetics; + + +//! Class that includes some bookeeping entries for a reaction or a global reaction defined on a surface +/*! + * Note that all indexes refer to a specific interfacial or homogeneous kinetics object. It does not + * refer to the Phase list indexes. + */ +class RxnMolChange +{ +public: + //! Main constructor for the class + /*! + * @param kinPtr Pointer to the kinetics base class + * @param irxn Specific reaction index. + */ + RxnMolChange(Cantera::Kinetics* kinPtr, int irxn); + + //! Destructor + ~RxnMolChange(); + + //! Constructor for the object if the object refers to a global reaction + /*! + * @param kinPtr Pointer to the kinetics base class + * @param egr Specific reaction index. + */ + RxnMolChange(Cantera::Kinetics* kinPtr, Cantera::ExtraGlobalRxn* egr); + + //! Vector of mole changes for each phase in the Kinetics object due to the current reaction + /*! + * This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient + * for all the species in a phase. + * The index is over all the phases listed in the Kinetics object. + */ + std::vector m_phaseMoleChange; + + std::vector m_phaseReactantMoles; + + std::vector m_phaseProductMoles; + + //! Vector of mass changes for each phase in the Kinetics object due to the current reaction + /*! + * This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient + * index multiplied by the molecular weight for all species in a phase. + * The index is over all of the phases listed in the Kinetics object. + */ + std::vector m_phaseMassChange; + + //! Vector of mass changes for each phase in the Kinetics object due to the current reaction + /*! + * This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient + * index multiplied by the charg for all species in a phase. + * The index is over all of the phases listed in the Kinetics object. + */ + std::vector m_phaseChargeChange; + + //! Vector of phase types in the reaction + /*! + * Collection of eosTypes for all phases in the kinetics object + * The index is over all of the phases listed in the Kinetics object. + */ + std::vector m_phaseTypes; + + //! Vector of phase dimensions for the reaction + /*! + * Collection of nDims for all phases in the kinetics object + * The index is over all of the phases listed in the Kinetics object. + */ + std::vector m_phaseDims; + + //! Number of phases in the kientics object + int m_nPhases; + + //! Shallow pointer pointing to the kinetics object + Cantera::Kinetics* m_kinBase; + + //! Reaction number within the kinetics object + /*! + * If this is neg 1, then this reaction refers to a global reaction + * specified by the m_egr pointer. + */ + int m_iRxn; + + //! Maximum change in charge of any phase due to this reaction + double m_ChargeTransferInRxn; + + //! Electrochemical beta parameter for the reaction + double m_beta; + + //! Pointer to the specification of the global reaction + /*! + * This is 0, if the class refers to a single reaction in the kinetics object + */ + Cantera::ExtraGlobalRxn* m_egr; +}; +} + +#endif diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h index dd9832408..7b33f83a3 100644 --- a/include/cantera/kinetics/StoichManager.h +++ b/include/cantera/kinetics/StoichManager.h @@ -8,6 +8,7 @@ #define CT_STOICH_MGR_H #include "cantera/base/stringUtils.h" +#include "cantera/base/ctexceptions.h" namespace Cantera { @@ -452,7 +453,8 @@ public: m_rxn(right.m_rxn), m_ic(right.m_ic), m_order(right.m_order), - m_stoich(right.m_stoich) { + m_stoich(right.m_stoich) + { } C_AnyN& operator=(const C_AnyN& right) { @@ -560,27 +562,44 @@ public: } private: + //! Length of the m_ic vector /*! - * This is the number of species which have non-zero entries in either the - * reaction order matrix or the stoichiometric order matrix for this reaction. + * This is the number of species which participate in the reaction order + * and stoichiometric coefficient vectors for the reactant or product description + * of the reaction. */ size_t m_n; //! ID of the reaction corresponding to this stoichiometric manager /*! - * This is used within the interface to select the + * This is used within the interface to select the array position to read and write to + * Normally this is associated with the reaction number in an array of quantities indexed + * by the reaction number, e.g., ROP[irxn]. */ size_t m_rxn; //! Vector of species which are involved with this stoichiometric manager calculations /*! - * This is an integer list of species which have non-zero entries in either the + * This is an integer list of species which participate in either the * reaction order matrix or the stoichiometric order matrix for this reaction, m_rxn. - * It's used as the index into the arrays m_order[] and m_stoich[]. */ std::vector m_ic; + + //! Reaction orders for the reaction + /*! + * This is either for the reactants or products. + * Length = m_n + * Species number, m_ic[n], has a reaction order of m_order[n]. + */ vector_fp m_order; + + //! Stoichiometric coefficients for the reaction, reactant or product side. + /*! + * This is either for the reactants or products. + * Length = m_n + * Species number m_ic[m], has a stoichiometric coefficient of m_stoich[n]. + */ vector_fp m_stoich; }; @@ -689,7 +708,7 @@ inline static void _writeMultiply(InputIter begin, InputIter end, * be the number of molecules on the product or reactant side of * reaction number i. * \f[ - * r_i = \sum_m^{M_i} s_{k_{m,i}} + * r_i = \sum_m^{M_i} s_{k_{m,i}} * \f] * To understand the operations performed by this class, let * \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on @@ -702,7 +721,7 @@ inline static void _writeMultiply(InputIter begin, InputIter end, * - \f$ S = S + N R\f$ (incrementSpecies) * - \f$ S = S - N R\f$ (decrementSpecies) * - \f$ R = R + N^T S \f$ (incrementReaction) - * - \f$ R = R - N^T S \f$ (deccrementReaction) + * - \f$ R = R - N^T S \f$ (decrementReaction) * * The actual implementation, however, does not compute these * quantities by matrix multiplication. A faster algorithm is used @@ -794,6 +813,12 @@ public: void add(size_t rxn, const std::vector& k, const vector_fp& order, const vector_fp& stoich) { m_n[rxn] = k.size(); + if (order.size() != k.size()) { + throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ"); + } + if (stoich.size() != k.size()) { + throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ"); + } bool frac = false; for (size_t n = 0; n < stoich.size(); n++) { if (stoich[n] != 1.0 || order[n] != 1.0) { @@ -900,11 +925,16 @@ private: std::vector m_c2_list; std::vector m_c3_list; std::vector m_cn_list; - /** - * Map with the Reaction Number as key and the Number of species - * as the value. + + //! Map with the Reaction Number as key and the Number of species + //! as the value. + /*! + * An example of this is given below: + * + * m_n[irxn] = nSpecies */ std::map m_n; + /** * Map with the Reaction Number as key and the placement in the * vector of reactions list( i.e., m_c1_list[]) as key diff --git a/src/kinetics/ElectrodeKinetics.cpp b/src/kinetics/ElectrodeKinetics.cpp new file mode 100644 index 000000000..37242ba40 --- /dev/null +++ b/src/kinetics/ElectrodeKinetics.cpp @@ -0,0 +1,74 @@ +/** + * @file ElectrodeKinetics.cpp + */ + +#include "cantera/kinetics/ElectrodeKinetics.h" + + +using namespace std; + +namespace Cantera +{ +//============================================================================================================================ +ElectrodeKinetics::ElectrodeKinetics(thermo_t* thermo) : + InterfaceKinetics(thermo), + metalPhaseRS_(npos), + electronPhaseRS_(npos), + solnPhaseRS_(npos), + kElectronRS_(npos) +{ + +} +//============================================================================================================================ +ElectrodeKinetics::~ElectrodeKinetics() +{ + +} +//============================================================================================================================ +ElectrodeKinetics::ElectrodeKinetics(const ElectrodeKinetics& right) : + InterfaceKinetics() + +{ + /* + * Call the assignment operator + */ + operator=(right); +} +//============================================================================================================================ +ElectrodeKinetics& ElectrodeKinetics::operator=(const ElectrodeKinetics& right) +{ + /* + * Check for self assignment. + */ + if (this == &right) { + return *this; + } + + InterfaceKinetics::operator=(right); + + metalPhaseRS_ = right.metalPhaseRS_; + electronPhaseRS_ = right.electronPhaseRS_; + solnPhaseRS_ = right.solnPhaseRS_; + kElectronRS_ = right.kElectronRS_; + + return *this; +} +//============================================================================================================================ +int ElectrodeKinetics::type() const +{ + return cInterfaceKinetics; +} +//============================================================================================================================ +Kinetics* ElectrodeKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const +{ + ElectrodeKinetics* iK = new ElectrodeKinetics(*this); + iK->assignShallowPointers(tpVector); + return iK; +} +//============================================================================================================================ + + + + +//================================================================================================================== +} diff --git a/src/kinetics/ExtraGlobalRxn.cpp b/src/kinetics/ExtraGlobalRxn.cpp new file mode 100644 index 000000000..0feccecca --- /dev/null +++ b/src/kinetics/ExtraGlobalRxn.cpp @@ -0,0 +1,394 @@ +/** + * @file example2.cpp + * + */ +/* + * $Id: ExtraGlobalRxn.cpp 571 2013-03-26 16:44:21Z hkmoffa $ + * + */ + +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +// Example 2 +// +// Read a mechanism, and print to the standard output stream a +// well-formatted Chemkin ELEMENT section. +// + +#include "cantera/kinetics/ExtraGlobalRxn.h" + + + +#include "cantera/numerics/DenseMatrix.h" + +// Kinetics includes +#include "cantera/kinetics.h" +#include "cantera/kinetics/InterfaceKinetics.h" +#include "cantera/thermo/SurfPhase.h" +#include "cantera/kinetics/KineticsFactory.h" + +#include +#include +#include +#include +#include + + + +using namespace std; +using namespace Cantera; + +namespace Cantera { + +//============================================================================================================ +static void erase_vd(std::vector& m_vec, int index) +{ + std::vector::iterator ipos; + ipos = m_vec.begin(); + ipos += index; + m_vec.erase(ipos); +} +//============================================================================================================ +static void erase_vi(std::vector& m_vec, int index) +{ + std::vector::iterator ipos; + ipos = m_vec.begin(); + ipos += index; + m_vec.erase(ipos); +} +//============================================================================================================ +//! add the species into the list of products or reactants +/*! + * Note this function gets called for both the product and reactant sides. However, it's only + * called for one side or another. + * + * @param kkinspec kinetic species index of the product + * @param + */ +static void addV(int kkinspec, double ps, std::vector& m_Products, + std::vector& m_ProductStoich) +{ + int nsize = m_Products.size(); + for (int i = 0; i < nsize; i++) { + if (m_Products[i] == kkinspec) { + m_ProductStoich[i] += ps; + return; + } + } + m_Products.push_back(kkinspec); + m_ProductStoich.push_back(ps); +} +//============================================================================================================ +ExtraGlobalRxn::ExtraGlobalRxn(Kinetics* k_ptr) : + m_ThisIsASurfaceRxn(false), + m_kinetics(k_ptr), + m_InterfaceKinetics(0), + m_nKinSpecies(0), + m_nReactants(0), + m_nProducts(0), + m_nNetSpecies(0), + m_nRxns(0), + m_SpecialSpecies(-1), + m_SpecialSpeciesProduct(true), + iphaseKin(0), + m_ok(false), + m_reversible(true) +{ + m_InterfaceKinetics = dynamic_cast(k_ptr); + if (m_InterfaceKinetics) { + m_ThisIsASurfaceRxn = true; + } + m_nRxns = m_kinetics->nReactions(); + m_ElemRxnVector.resize(m_nRxns,0.0); + m_nKinSpecies = m_kinetics->nTotalSpecies(); +} +//============================================================================================================ +ExtraGlobalRxn::~ExtraGlobalRxn() +{ +} +//============================================================================================================ +void ExtraGlobalRxn::setupElemRxnVector(double* RxnVector, + int specialSpecies) +{ + int i; + int kkinspec; + for (size_t i = 0; i < (size_t) m_nRxns; i++) { + m_ElemRxnVector[i] = RxnVector[i]; + } + for (size_t i = 0; i < (size_t) m_nRxns; i++) { + if (RxnVector[i] > 0.0) { + for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) { + double ps = m_kinetics->productStoichCoeff(kkinspec, i); + if (ps > 0.0) { + addV(kkinspec, RxnVector[i]* ps, m_Products, m_ProductStoich); + addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich); + } + double rs = m_kinetics->reactantStoichCoeff(kkinspec, i); + if (rs > 0.0) { + addV(kkinspec, RxnVector[i] * rs, m_Reactants, m_ReactantStoich); + addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich); + } + } + } else if (RxnVector[i] < 0.0) { + for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) { + double ps = m_kinetics->productStoichCoeff(kkinspec, i); + if (ps > 0.0) { + addV(kkinspec,- RxnVector[i]* ps, m_Reactants, m_ReactantStoich); + addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich); + } + double rs = m_kinetics->reactantStoichCoeff(kkinspec, i); + if (rs > 0.0) { + addV(kkinspec, -RxnVector[i] * rs, m_Products, m_ProductStoich); + addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich); + } + } + } + } +Recheck: + for (i = 0; i < static_cast(m_Products.size()); i++) { + if (m_ProductStoich[i] == 0.0) { + erase_vi(m_Products, i); + erase_vd(m_ProductStoich, i); + goto Recheck ; + } + } + for (i = 0; i < static_cast(m_Reactants.size()); i++) { + if (m_ReactantStoich[i] == 0.0) { + erase_vi(m_Reactants, i); + erase_vd(m_ReactantStoich, i); + goto Recheck ; + } + } + for (i = 0; i < static_cast(m_NetSpecies.size()); i++) { + if (m_netStoich[i] == 0.0) { + erase_vi(m_NetSpecies, i); + erase_vd(m_netStoich, i); + goto Recheck ; + } + } + + for (i = 0; i < static_cast(m_Products.size()); i++) { + int ik = m_Products[i]; + for (int j = 0; j < static_cast(m_Reactants.size()); j++) { + int jk = m_Reactants[j]; + if (ik == jk) { + if (m_ProductStoich[i] == m_ReactantStoich[j]) { + erase_vi(m_Products, i); + erase_vd(m_ProductStoich, i); + erase_vi(m_Reactants, j); + erase_vd(m_ReactantStoich, j); + } else if (m_ProductStoich[i] > m_ReactantStoich[j]) { + m_ProductStoich[i] -= m_ReactantStoich[j]; + erase_vi(m_Reactants, j); + erase_vd(m_ReactantStoich, j); + } else { + m_ReactantStoich[j] -= m_ProductStoich[i]; + erase_vi(m_Products, i); + erase_vd(m_ProductStoich, i); + } + // We just screwed up the indexing -> restart. + goto Recheck ; + } + + } + } + m_nProducts = m_Products.size(); + m_nReactants = m_Reactants.size(); + m_nNetSpecies = m_NetSpecies.size(); + + /* + * Section to assign the special species + */ + m_SpecialSpecies = specialSpecies; + if (specialSpecies == -1) { + m_SpecialSpecies = m_Products[0]; + } + bool ifound = false; + for (i = 0; i < (int) m_NetSpecies.size(); i++) { + int ik = m_NetSpecies[i]; + if (ik == m_SpecialSpecies) { + if (m_netStoich[i] > 0.0) { + m_SpecialSpeciesProduct = true; + } else { + m_SpecialSpeciesProduct = false; + } + m_SS_index = i; + ifound = true; + break; + } + } + if (!ifound) { + throw CanteraError(":setupElemRxnVector", + "Special species not a reactant or product: " + + int2str(m_SpecialSpecies)); + } + + m_ok = true; +} +//============================================================================================================ +std::string ExtraGlobalRxn::reactionString() +{ + string rs; + int k, istoich; + for (k = 0; k < m_nReactants; k++) { + int kkinspecies = m_Reactants[k]; + double stoich = m_ReactantStoich[k]; + if (stoich != 1.0) { + istoich = (int) stoich; + if (fabs((double)istoich - stoich) < 0.00001) { + rs += int2str(istoich) + " "; + } else { + rs += fp2str(stoich) + " "; + } + } + string sName = m_kinetics->kineticsSpeciesName(kkinspecies); + rs += sName; + if (k < (m_nReactants-1)) { + rs += " + "; + } + } + rs += " = "; + for (k = 0; k < m_nProducts; k++) { + int kkinspecies = m_Products[k]; + double stoich = m_ProductStoich[k]; + if (stoich != 1.0) { + istoich = (int) stoich; + if (fabs((double)istoich - stoich) < 0.00001) { + rs += int2str(istoich) + " "; + } else { + rs += fp2str(stoich) + " "; + } + } + string sName = m_kinetics->kineticsSpeciesName(kkinspecies); + rs += sName; + if (k < (m_nProducts-1)) { + rs += " + "; + } + } + return rs; +} +//============================================================================================================ + +std::vector& ExtraGlobalRxn::reactants() +{ + return m_Reactants; +} +//============================================================================================================ + +std::vector& ExtraGlobalRxn::products() +{ + return m_Products; +} + +//============================================================================================================ +bool ExtraGlobalRxn::isReversible() +{ + return m_reversible; +} +//============================================================================================================ + +double ExtraGlobalRxn::reactantStoichCoeff(int kKin) +{ + for (int k = 0; k < m_nReactants; k++) { + int kkinspec = m_Reactants[k]; + if (kkinspec == kKin) { + return m_ReactantStoich[k]; + } + } + return 0.0; +} +//============================================================================================================ +double ExtraGlobalRxn::productStoichCoeff(int kKin) +{ + for (int k = 0; k < m_nProducts; k++) { + int kkinspec = m_Products[k]; + if (kkinspec == kKin) { + return m_ProductStoich[k]; + } + } + return 0.0; +} +//============================================================================================================ + +double ExtraGlobalRxn::deltaSpecValue(double* speciesVectorProperty) +{ + int k; + double val = 0; + for (k = 0; k < m_nNetSpecies; k++) { + int kkinspec = m_NetSpecies[k]; + val += speciesVectorProperty[kkinspec] * m_netStoich[k]; + } + return val; +} +//============================================================================================================ + +double ExtraGlobalRxn::deltaRxnVecValue(double* rxnVectorProperty) +{ + double val = 0; + for (int i = 0; i < m_nRxns; i++) { + val += m_ElemRxnVector[i] * rxnVectorProperty[i]; + } + return val; +} +//============================================================================================================ +double ExtraGlobalRxn::ROPValue(double* ROPElemKinVector) +{ + double val = 0.0; + for (int i = 0; i < m_nRxns; i++) { + double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i); + if (m_ElemRxnVector[i] > 0.0) { + val += ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } else { + val -= ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } + } + if (!m_SpecialSpeciesProduct) { + val = -val; + } + return val; +} +//============================================================================================================ +double ExtraGlobalRxn::FwdROPValue(double* FwdROPElemKinVector, + double* RevROPElemKinVector) +{ + double val = 0.0; + for (int i = 0; i < m_nRxns; i++) { + double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i); + if (m_ElemRxnVector[i] > 0.0) { + val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } + if (m_ElemRxnVector[i] < 0.0) { + val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } + } + if (!m_SpecialSpeciesProduct) { + val = -val; + } + return val; +} +//============================================================================================================ +double ExtraGlobalRxn::RevROPValue(double* FwdROPElemKinVector, + double* RevROPElemKinVector) +{ + double val = 0.0; + for (int i = 0; i < m_nRxns; i++) { + double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i)- m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i); + if (m_ElemRxnVector[i] > 0.0) { + val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } + if (m_ElemRxnVector[i] < 0.0) { + val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i]; + } + } + if (!m_SpecialSpeciesProduct) { + val = -val; + } + return val; +} +//============================================================================================================ +} + diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 4d81968e1..44e86f2f7 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -31,6 +31,7 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) : m_ctrxn_ecdf(0), m_StandardConc(0), m_deltaG0(0), + m_deltaG(0), m_ProdStanConcReac(0), m_logp0(0.0), m_logc0(0.0), @@ -71,6 +72,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : m_ctrxn_ecdf(0), m_StandardConc(0), m_deltaG0(0), + m_deltaG(0), m_ProdStanConcReac(0), m_logp0(0.0), m_logc0(0.0), @@ -121,6 +123,7 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) m_conc = right.m_conc; m_actConc = right.m_actConc; m_mu0 = right.m_mu0; + m_mu = right.m_mu; m_mu0_Kc = right.m_mu0_Kc; m_phi = right.m_phi; m_pot = right.m_pot; @@ -134,6 +137,7 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) m_ctrxn_ecdf = right.m_ctrxn_ecdf; m_StandardConc = right.m_StandardConc; m_deltaG0 = right.m_deltaG0; + m_deltaG = right.m_deltaG; m_ProdStanConcReac = right.m_ProdStanConcReac; m_logp0 = right.m_logp0; m_logc0 = right.m_logc0; @@ -347,13 +351,13 @@ void InterfaceKinetics::getRevRatesOfProgress(doublereal* revROP) updateROP(); std::copy(m_ropr.begin(), m_ropr.end(), revROP); } - +//=========================================================================================================== void InterfaceKinetics::getNetRatesOfProgress(doublereal* netROP) { updateROP(); std::copy(m_ropnet.begin(), m_ropnet.end(), netROP); } - +//=========================================================================================================== void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) { updateMu0(); @@ -367,8 +371,9 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) kc[i] = exp(-kc[i]*rrt); } } - -/** values needed to convert from exchange current density to surface reaction rate. +//=========================================================================================================== +/* + * values needed to convert from exchange current density to surface reaction rate. */ void InterfaceKinetics::updateExchangeCurrentQuantities() { @@ -391,26 +396,28 @@ void InterfaceKinetics::updateExchangeCurrentQuantities() m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_deltaG0)); + // + // Calculate the product of the standard concentrations of the reactants + // for (size_t i = 0; i < m_ii; i++) { m_ProdStanConcReac[i] = 1.0; } - m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac)); } - +//=========================================================================================================== void InterfaceKinetics::getCreationRates(doublereal* cdot) { updateROP(); m_rxnstoich.getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot); } - +//=========================================================================================================== void InterfaceKinetics::getDestructionRates(doublereal* ddot) { updateROP(); m_rxnstoich.getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot); } - +//=========================================================================================================== void InterfaceKinetics::getNetProductionRates(doublereal* net) { updateROP(); @@ -566,34 +573,87 @@ void InterfaceKinetics::updateROP() _update_rates_T(); // get updated activities (rates updated below) _update_rates_C(); + + double TT = m_surf->temperature(); + double rtdf = GasConstant * TT / Faraday; if (m_ROP_ok) { return; } - - // copy rate coefficients into ropf + // + // Copy the reaction rate coefficients, m_rfn, into m_ropf + // copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); - - // multiply by perturbation factor + // + // Multiply by the perturbation factor + // multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin()); - - // copy the forward rates to the reverse rates + // + // Copy the forward rate constants to the reverse rate constants + // copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin()); - // for reverse rates computed from thermochemistry, multiply + + + // + // For reverse rates computed from thermochemistry, multiply // the forward rates copied into m_ropr by the reciprocals of // the equilibrium constants + // multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); // // multiply ropf by the actyivity concentration reaction orders to obtain // the forward rates of progress. // m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf)); + // + // For reversible reactions, multiply ropr by the activity concentration products + // + m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc), DATA_PTR(m_ropr)); + // + // Fix up these calculations for cases where the above formalism doesn't hold + // + double OCV = 0.0; + for (size_t jrxn = 0; jrxn != m_ii; ++jrxn) { + int reactionType = reactionTypes_[jrxn]; + if (reactionType == BUTLERVOLMER_RXN) { + // + // OK, the reaction rate constant contains the current density rate constant calculation + // the rxnstoich calculation contained the dependence of the current density on the activity concentrations + // We finish up with the ROP calculation + // + // Calculate the overpotential of the reaction + // + // double nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_]; + double nStoichElectrons=1; + //*nStoich = nStoichElectrons; + + + + getDeltaGibbs(0); + + if (nStoichElectrons != 0.0) { + OCV = m_deltaG[jrxn]/Faraday/ nStoichElectrons; + } + + /* + + double exp1 = nu * nStoich * beta / rtdf + double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (GasConstant * temp); + double val = io * (exp(exp1) - exp(exp2)); + + doublereal BVterm = exp(exp1 ) - exp(exp2); + m_ropnet[j] = m_ropf[j] * BVterm + m_ropf[j] = + // + m_ropr[j] = m_ropnet[j] - m_ropf[j]; + */ + } + + + } + - // for reversible reactions, multiply ropr by the activity concentration - // products - m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc), - DATA_PTR(m_ropr)); for (size_t j = 0; j != m_ii; ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -670,13 +730,18 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) * kinetics mechanism */ for (size_t n = 0; n < nPhases(); n++) { - thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]); + m_thermo[n]->getChemPotentials(DATA_PTR(m_mu) + m_start[n]); + } + // + // Use the stoichiometric manager to find deltaG for each + // reaction. + // + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu), DATA_PTR(m_deltaG)); + if (deltaG != 0 && (DATA_PTR(m_deltaG) != deltaG)) { + for (size_t j = 0; j < m_ii; ++j) { + deltaG[j] = m_deltaG[j]; + } } - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); } //================================================================================================================== void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) @@ -709,7 +774,7 @@ void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); } - +//=========================================================================================================== void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS) { /* @@ -725,7 +790,7 @@ void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } - +//=========================================================================================================== void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS) { /* @@ -743,7 +808,7 @@ void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), deltaGSS); } - +//=========================================================================================================== void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) { /* @@ -755,7 +820,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) for (size_t n = 0; n < nPhases(); n++) { thermo(n).getEnthalpy_RT(DATA_PTR(m_grt) + m_start[n]); } - doublereal RT = thermo().temperature() * GasConstant; + doublereal RT = thermo(0).temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { m_grt[k] *= RT; } @@ -765,7 +830,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); } - +//=========================================================================================================== void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS) { /* @@ -859,6 +924,12 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) rp.push_back(rdata.cov[m]); } + // + // Find out the reaction type + // + int reactionType = rdata.reactionType; + reactionTypes_.push_back(reactionType); + /* * Temporarily change the reaction rate coefficient type to surface arrhenius. * This is what is expected. We'll handle exchange current types below by hand. @@ -919,7 +990,8 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) // // Find out the reaction type // - int reactionType = rdata.reactionType; + int reactionType = rdata.reactionType; + reactionTypes_.push_back(reactionType); /* * Temporarily change the reaction rate coefficient type to surface arrhenius. @@ -945,10 +1017,15 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) // store activation energy m_E.push_back(rdata.rateCoeffParameters[2]); + // + // Add the reaction into the list of electrochemical extras + // if (rdata.beta > 0.0) { m_has_electrochem_rxns = true; m_beta.push_back(rdata.beta); + // Push back the id of the reaction m_ctrxn.push_back(m_ii); + // set the default to be the normal forward / reverse calculation method m_ctrxn_BVform.push_back(0); if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) { m_has_exchange_current_density_formulation = true; @@ -1083,6 +1160,7 @@ void InterfaceKinetics::init() m_actConc.resize(m_kk); m_conc.resize(m_kk); m_mu0.resize(m_kk); + m_mu.resize(m_kk); m_mu0_Kc.resize(m_kk); m_grt.resize(m_kk); m_pot.resize(m_kk, 0.0); @@ -1105,6 +1183,8 @@ void InterfaceKinetics::finalize() m_StandardConc.resize(m_kk, 0.0); m_deltaG0.resize(safe_reaction_size, 0.0); + m_deltaG.resize(safe_reaction_size, 0.0); + m_ProdStanConcReac.resize(safe_reaction_size, 0.0); if (m_thermo.size() != m_phaseExists.size()) { @@ -1219,6 +1299,20 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl } } //================================================================================================================== +// Write values into m_index +/* + * @param rxnNumber reaction number + * @param type reaction type + * @param loc location location in the reaction rate coefficient calculator + */ +void InterfaceKinetics::registerReaction(size_t rxnNumber, int type, size_t loc) +{ + // + // type and loc is storred as a pair of values. + // + m_index[rxnNumber] = std::pair(type, loc); +} +//================================================================================================================== void EdgeKinetics::finalize() { deltaElectricEnergy_.resize(std::max(m_ii, 1)); @@ -1243,5 +1337,6 @@ void EdgeKinetics::finalize() m_finalized = true; } + //================================================================================================================== } diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index f52cee605..e70b3fbfb 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -1,7 +1,6 @@ /** * @file Kinetics.cpp Declarations for the base class for kinetics managers - * (see \ref kineticsmgr and class \link Cantera::Kinetics - * Kinetics\endlink). + * (see \ref kineticsmgr and class \link Cantera::Kinetics Kinetics \endlink). * * Kinetics managers calculate rates of progress of species due to * homogeneous or heterogeneous kinetics. diff --git a/src/kinetics/ReactionStoichMgr.cpp b/src/kinetics/ReactionStoichMgr.cpp index b99a30702..695e147b9 100644 --- a/src/kinetics/ReactionStoichMgr.cpp +++ b/src/kinetics/ReactionStoichMgr.cpp @@ -45,7 +45,7 @@ ReactionStoichMgr& ReactionStoichMgr::operator=(const ReactionStoichMgr& right) } return *this; } - +//========================================================================================================= void ReactionStoichMgr::add(size_t rxn, const std::vector& reactants, const std::vector& products, bool reversible) @@ -59,12 +59,65 @@ void ReactionStoichMgr::add(size_t rxn, const std::vector& reactants, m_irrevproducts.add(rxn, products); } } - +//========================================================================================================= +// Add the reaction into the stoichiometric manager void ReactionStoichMgr::add(size_t rxn, const ReactionData& r) { - + size_t k; std::vector rk; doublereal frac; + doublereal oo, os, of; + bool doGlobal = false; + std::vector extReactants = r.reactants; + vector_fp extRStoich = r.rstoich; + vector_fp extROrder = r.rorder; + + // + // If we have a complete global reaction then we need to do something more complete + // than the previous treatment. Basically we will use the reactant manager to calculate the + // global forward reaction rate of progress. + // + if (r.forwardFullOrder_.size() > 0) { + // + // Trigger a treatment where the order of the reaction and the stoichiometry + // are treated as different. + // + doGlobal = true; + size_t nsp = r.forwardFullOrder_.size(); + // + // Set up a signal vector to indicate whether the species has been added into + // the input vectors for the stoich manager + // + vector_int kHandled(nsp, 0); + // + // Loop over the reactants which are also nonzero stoichioemtric entries + // making sure the forwardFullOrder_ entries take precedence over rorder entries + // + for (size_t kk = 0; kk < r.reactants.size(); kk++) { + k = r.reactants[kk]; + os = r.rstoich[kk]; + oo = r.rorder[kk]; + of = r.forwardFullOrder_[k]; + if (of != oo) { + extROrder[kk] = of; + } + kHandled[k] = 1; + } + for (k = 0; k < nsp; k++) { + of = r.forwardFullOrder_[k]; + if (of != 0.0) { + if (kHandled[k] == 0) { + // + // Add extra entries to reactant inputs. Set their reactant stoichiometric entries to zero. + // + extReactants.push_back(k); + extROrder.push_back(of); + extRStoich.push_back(0.0); + } + } + } + } + bool isfrac = false; for (size_t n = 0; n < r.reactants.size(); n++) { size_t ns = size_t(r.rstoich[n]); @@ -77,12 +130,22 @@ void ReactionStoichMgr::add(size_t rxn, const ReactionData& r) } } - // if the reaction has fractional stoichiometric coefficients - // or specified reaction orders, then add it in a general reaction - if (isfrac || r.global || rk.size() > 3) { - m_reactants.add(rxn, r.reactants, r.rorder, r.rstoich); + // + // If the reaction is non-mass action add it in in a general way + // Reactants get extra terms for the forward reaction rate of progress + // that may have zero stoichiometries. + // + if (doGlobal) { + m_reactants.add(rxn, extReactants, extROrder, extRStoich); } else { - m_reactants.add(rxn, rk); + // + // this is confusing. The only issue should be whether rorder is different than rstoich! + // + if (isfrac || r.global || rk.size() > 3) { + m_reactants.add(rxn, r.reactants, r.rorder, r.rstoich); + } else { + m_reactants.add(rxn, rk); + } } std::vector pk; @@ -99,22 +162,26 @@ void ReactionStoichMgr::add(size_t rxn, const ReactionData& r) } if (r.reversible) { - if (isfrac && !r.isReversibleWithFrac) { - throw CanteraError("ReactionStoichMgr::add", - "Fractional product stoichiometric coefficients only allowed " - "\nfor irreversible reactions and most reversible reactions"); - } + // + // this is confusing. The only issue should be whether porder is different than pstoich! + // if (pk.size() > 3 || r.isReversibleWithFrac) { m_revproducts.add(rxn, r.products, r.porder, r.pstoich); } else { m_revproducts.add(rxn, pk); } - } else if (isfrac || pk.size() > 3) { - m_irrevproducts.add(rxn, r.products, r.porder, r.pstoich); } else { - m_irrevproducts.add(rxn, pk); + // + // this is confusing. The only issue should be whether porder is different than pstoich! + // + if (isfrac || pk.size() > 3) { + m_irrevproducts.add(rxn, r.products, r.porder, r.pstoich); + } else { + m_irrevproducts.add(rxn, pk); + } } } +//========================================================================================================= void ReactionStoichMgr::getCreationRates(size_t nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c) diff --git a/src/kinetics/RxnMolChange.cpp b/src/kinetics/RxnMolChange.cpp new file mode 100644 index 000000000..ced443165 --- /dev/null +++ b/src/kinetics/RxnMolChange.cpp @@ -0,0 +1,164 @@ +/** + * @file example2.cpp + * + * $Id: RxnMolChange.cpp 571 2013-03-26 16:44:21Z hkmoffa $ + * + */ + +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#include "cantera/kinetics/RxnMolChange.h" + + +#include "cantera/thermo.h" +#include "cantera/kinetics.h" +#include "cantera/kinetics/InterfaceKinetics.h" + +#include "cantera/kinetics/ExtraGlobalRxn.h" + +#include +#include + +using namespace Cantera; +using namespace std; + +namespace Cantera { + +RxnMolChange::RxnMolChange(Cantera::Kinetics* kinPtr, int irxn) : + m_nPhases(0), + m_kinBase(kinPtr), + m_iRxn(irxn), + m_ChargeTransferInRxn(0.0), + m_beta(0.0), + m_egr(0) +{ + int nReac = kinPtr->nReactions(); + int iph; + AssertTrace(irxn >= 0); + AssertTrace(irxn < nReac); + + m_nPhases = kinPtr->nPhases(); + + m_phaseMoleChange.resize(m_nPhases, 0.0); + m_phaseReactantMoles.resize(m_nPhases, 0.0); + m_phaseProductMoles.resize(m_nPhases, 0.0); + m_phaseMassChange.resize(m_nPhases, 0.0); + m_phaseChargeChange.resize(m_nPhases, 0.0); + m_phaseTypes.resize(m_nPhases, 0); + m_phaseDims.resize(m_nPhases, 0); + + int m_kk = kinPtr->nTotalSpecies(); + + for (int kKin = 0; kKin < m_kk; kKin++) { + iph = m_kinBase->speciesPhaseIndex(kKin); + Cantera::ThermoPhase& tpRef = m_kinBase->thermo(iph); + int kLoc = kKin - m_kinBase->kineticsSpeciesIndex(0, iph); + double rsc = m_kinBase->reactantStoichCoeff(kKin, irxn); + double psc = m_kinBase->productStoichCoeff(kKin, irxn); + double nsc = psc - rsc; + m_phaseMoleChange[iph] += (nsc); + m_phaseReactantMoles[iph] += rsc; + m_phaseProductMoles[iph] += psc; + double mw = tpRef.molecularWeight(kLoc); + m_phaseMassChange[iph] += (nsc) * mw; + double chg = tpRef.charge(kLoc); + m_phaseChargeChange[iph] += nsc * chg; + } + + for (iph = 0; iph < m_nPhases; iph++) { + Cantera::ThermoPhase& tpRef = m_kinBase->thermo(iph); + m_phaseDims[iph] = tpRef.nDim(); + m_phaseTypes[iph] = tpRef.eosType(); + if (m_phaseChargeChange[iph] != 0.0) { + double tmp = fabs(m_phaseChargeChange[iph]); + if (tmp > m_ChargeTransferInRxn) { + m_ChargeTransferInRxn = tmp; + } + } + } + + if (m_ChargeTransferInRxn) { + Cantera::InterfaceKinetics* iK = dynamic_cast(kinPtr); + if (iK) { + m_beta = iK->electrochem_beta(irxn); + } else { + throw Cantera::CanteraError("RxnMolChange", "unknown condition on charge"); + } + } + +} + +RxnMolChange::RxnMolChange(Cantera::Kinetics* kinPtr, Cantera::ExtraGlobalRxn* egr) : + m_nPhases(0), + m_kinBase(kinPtr), + m_iRxn(-1), + m_ChargeTransferInRxn(0.0), + m_beta(0.0), + m_egr(egr) +{ + int iph; + AssertTrace(egr != 0); + + m_nPhases = kinPtr->nPhases(); + + m_phaseMoleChange.resize(m_nPhases, 0.0); + m_phaseReactantMoles.resize(m_nPhases, 0.0); + m_phaseProductMoles.resize(m_nPhases, 0.0); + m_phaseMassChange.resize(m_nPhases, 0.0); + m_phaseChargeChange.resize(m_nPhases, 0.0); + m_phaseTypes.resize(m_nPhases, 0); + m_phaseDims.resize(m_nPhases, 0); + + int m_kk = kinPtr->nTotalSpecies(); + + for (int kKin = 0; kKin < m_kk; kKin++) { + iph = m_kinBase->speciesPhaseIndex(kKin); + ThermoPhase& tpRef = m_kinBase->thermo(iph); + int kLoc = kKin - m_kinBase->kineticsSpeciesIndex(0, iph); + double rsc = egr->reactantStoichCoeff(kKin); + double psc = egr->productStoichCoeff(kKin); + double nsc = psc - rsc; + m_phaseMoleChange[iph] += (nsc); + m_phaseReactantMoles[iph] += rsc; + m_phaseProductMoles[iph] += psc; + double mw = tpRef.molecularWeight(kLoc); + m_phaseMassChange[iph] += (nsc) * mw; + double chg = tpRef.charge(kLoc); + m_phaseChargeChange[iph] += nsc * chg; + } + + for (iph = 0; iph < m_nPhases; iph++) { + ThermoPhase& tpRef = m_kinBase->thermo(iph); + m_phaseDims[iph] = tpRef.nDim(); + m_phaseTypes[iph] = tpRef.eosType(); + if (m_phaseChargeChange[iph] != 0.0) { + double tmp = fabs(m_phaseChargeChange[iph]); + if (tmp > m_ChargeTransferInRxn) { + m_ChargeTransferInRxn = tmp; + } + } + } + + if (m_ChargeTransferInRxn) { + InterfaceKinetics* iK = dynamic_cast(kinPtr); + if (iK) { + m_beta = 0.0; + } else { + throw CanteraError("RxnMolChange", "unknown condition on charge"); + } + } + +} + + +RxnMolChange::~RxnMolChange() +{ + +} + +} + diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 389c55bc6..d8e74ef0c 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -82,7 +82,7 @@ public: } } }; - +//======================================================================================================= void checkRxnElementBalance(Kinetics& kin, const ReactionData& rdata, doublereal errorTolerance) {