From e7236df51f6e6a5e55dbea761fa82579b0aaedad Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 22 Oct 2014 21:48:09 +0000 Subject: [PATCH] Formatting cleanup of InterfaceKinetics --- include/cantera/kinetics/InterfaceKinetics.h | 43 +- src/kinetics/InterfaceKinetics.cpp | 648 ++++++++----------- 2 files changed, 295 insertions(+), 396 deletions(-) diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 3c58f03ff..1cb4b1c03 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -19,12 +19,12 @@ namespace Cantera { -// forward references +// forward declarations class SurfPhase; class ImplicitSurfChem; class RxnMolChange; -//! forward orders +//! forward orders class RxnOrders { public: @@ -36,18 +36,18 @@ class RxnOrders { ~RxnOrders(); RxnOrders& operator=(const RxnOrders &right); - + //! Fill in the structure with the array. /*! * @param[in] Size of length kinetic species. The entries the values of the orders */ - int fill(const std::vector& fullForwardOrders); + int fill(const vector_fp& fullForwardOrders); //! ID's of the kinetic species std::vector kinSpeciesIDs_; //! Orders of the kinetic species - std::vector kinSpeciesOrders_; + vector_fp kinSpeciesOrders_; }; //! A kinetics manager for heterogeneous reaction mechanisms. The @@ -187,9 +187,9 @@ public: /*! * @param[in] Reaction index * - * @return Returns the reaction type of the reaction. + * @return Returns the reaction type of the reaction. */ - virtual int reactionType(size_t irxn) const; + virtual int reactionType(size_t irxn) const; virtual void getActivityConcentrations(doublereal* const conc); @@ -336,16 +336,14 @@ public: //! Add a single elementary reaction to the list of reactions for the object /*! - * @param rdata + * @param rdata */ void addElementaryReaction(ReactionData& rdata); - void addGlobalReaction(ReactionData& r); void installReagents(const ReactionData& r); - //! Update the equilibrium constants and stored electrochemical potentials //! in molar units for all reversible reactions and for all species. /*! @@ -390,7 +388,7 @@ public: * reaction type is a butler-volmer form, convert it to exchange current density * form (amps/m2). * - * @param kfwd Vector of forward reaction rate constants, given in either + * @param kfwd Vector of forward reaction rate constants, given in either * normal form or in exchange current density form. */ void convertExchangeCurrentDensityFormulation(doublereal* const kfwd); @@ -447,15 +445,12 @@ public: */ int phaseStability(const size_t iphase) const; - - virtual void determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders); + virtual void determineFwdOrdersBV(ReactionData& rdata, vector_fp& fwdFullorders); 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 @@ -523,7 +518,7 @@ protected: //! Vector of reactionType for the reactions defined within this object /*! * Length = number of reactions, m_ii - * contains the type of reaction. + * contains the type of reaction. */ vector_int reactionType_; @@ -595,7 +590,7 @@ protected: * potentials + RTln(Cs) for all of the species in the kinetics object * * In order to get the units correct for the concentration equilibrium - * constant, each species needs to have an + * constant, each species needs to have an * RT ln(Cs) added to its contribution to the equilibrium constant * Cs is the standard concentration for the species. Frequently, for * solid species, Cs is equal to 1. However, for gases Cs is P/RT. @@ -670,12 +665,12 @@ protected: */ std::vector m_ctrxn; - //! Vector of Reactions which follow the butler volmer methodology for specifying the + //! Vector of Reactions which follow the butler volmer methodology for specifying the //! exchange current density first. Then, the other forms are specified based on this form. /*! * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] - * - * m_ctrxn_BVform[i] = 0; This means that the irxn reaction is calculated via the standard forward + * + * m_ctrxn_BVform[i] = 0; This means that the irxn reaction is calculated via the standard forward * and reverse reaction rates * m_ctrxn_BVform[i] = 1; This means that the irxn reaction is calculated via the BV format * directly. @@ -684,7 +679,7 @@ protected: */ std::vector m_ctrxn_BVform; - //! Vector of booleans indicating whether the charge transfer reaction rate constant + //! Vector of booleans indicating whether the charge transfer reaction rate constant //! is described by an exchange current density rate constant expression /*! * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] @@ -696,7 +691,7 @@ protected: */ vector_int m_ctrxn_ecdf; - //! Vector of booleans indicating whether the charge transfer reaction rate constant + //! Vector of booleans indicating whether the charge transfer reaction rate constant //! is described by an exchange current density rate constant expression /*! * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] @@ -713,7 +708,7 @@ protected: */ std::vector m_ctrxn_FwdOrdersList_; - std::vector m_ctrxn_resistivity_; + vector_fp m_ctrxn_resistivity_; //! Vector of standard concentrations /*! @@ -811,7 +806,7 @@ protected: /*! * 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. - * For example, a reaction can't go in the forwards direction if a + * For example, a reaction can't go in the forwards direction if a * phase in which a reactant is present doesn't exist. Because InterfaceKinetics * deals with intrinsic quantities only normally, nowhere else is this extrinsic * concept introduced except here. diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 8eb003c16..53256fcc2 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -17,7 +17,7 @@ using namespace std; namespace Cantera { -//============================================================================================================================ + InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) : Kinetics(), m_redo_rates(false), @@ -53,7 +53,7 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) : addPhase(*thermo); } } -//============================================================================================================================ + InterfaceKinetics::~InterfaceKinetics() { delete m_integrator; @@ -68,7 +68,7 @@ InterfaceKinetics::~InterfaceKinetics() delete m_ctrxn_FwdOrdersList_[i]; } } -//============================================================================================================================ + InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : Kinetics(), m_redo_rates(false), @@ -105,7 +105,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : */ operator=(right); } -//============================================================================================================================ + InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) { /* @@ -185,8 +185,8 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) } m_ctrxn_ROPOrdersList_ = right.m_ctrxn_ROPOrdersList_; for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) { - RxnOrders* ro = right.m_ctrxn_ROPOrdersList_[i]; - m_ctrxn_ROPOrdersList_[i] = new RxnOrders(*ro); + RxnOrders* ro = right.m_ctrxn_ROPOrdersList_[i]; + m_ctrxn_ROPOrdersList_[i] = new RxnOrders(*ro); } for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) { @@ -194,59 +194,52 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) } m_ctrxn_FwdOrdersList_ = right.m_ctrxn_FwdOrdersList_; for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) { - RxnOrders* ro = right.m_ctrxn_FwdOrdersList_[i]; - m_ctrxn_FwdOrdersList_[i] = new RxnOrders(*ro); + RxnOrders* ro = right.m_ctrxn_FwdOrdersList_[i]; + m_ctrxn_FwdOrdersList_[i] = new RxnOrders(*ro); } - return *this; } -//============================================================================================================================ + int InterfaceKinetics::type() const { return cInterfaceKinetics; } -//============================================================================================================================ + Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const { InterfaceKinetics* iK = new InterfaceKinetics(*this); iK->assignShallowPointers(tpVector); return iK; } -//============================================================================================================================ + void InterfaceKinetics::setElectricPotential(int n, doublereal V) { thermo(n).setElectricPotential(V); m_redo_rates = true; } -//============================================================================================================================ + void InterfaceKinetics::_update_rates_T() { - // // First task is update the electrical potentials from the Phases - // _update_rates_phi(); if (m_has_coverage_dependence) { m_surf->getCoverages(DATA_PTR(m_actConc)); m_rates.update_C(DATA_PTR(m_actConc)); m_redo_rates = true; } - // + // Go find the temperature from the surface - // doublereal T = thermo(surfacePhaseIndex()).temperature(); m_redo_rates = true; if (T != m_temp || m_redo_rates) { - m_logtemp = log(T); - // - // Calculate the forward rate constant by calling m_rates and store it in m_rfn[] - // + + // Calculate the forward rate constant by calling m_rates and store it in m_rfn[] m_rates.update(T, m_logtemp, DATA_PTR(m_rfn)); - // - // If we need to do conversions between exchange current density formulation and regular formulation - // (either way) do it here. - // + + // If we need to do conversions between exchange current density formulation and regular formulation + // (either way) do it here. if (m_has_exchange_current_density_formulation) { convertExchangeCurrentDensityFormulation(DATA_PTR(m_rfn)); } @@ -259,12 +252,10 @@ void InterfaceKinetics::_update_rates_T() m_redo_rates = false; } } -//============================================================================================================================ + void InterfaceKinetics::_update_rates_phi() { - // // Store electric potentials for each phase in the array m_phi[]. - // for (size_t n = 0; n < nPhases(); n++) { if (thermo(n).electricPotential() != m_phi[n]) { m_phi[n] = thermo(n).electricPotential(); @@ -272,14 +263,12 @@ void InterfaceKinetics::_update_rates_phi() } } } -//============================================================================================================================ -// + // Updates the internal variables m_actConc and m_conc -// void InterfaceKinetics::_update_rates_C() { for (size_t n = 0; n < nPhases(); n++) { - const ThermoPhase* tp = m_thermo[n]; + const ThermoPhase* tp = m_thermo[n]; /* * We call the getActivityConcentrations function of each * ThermoPhase class that makes up this kinetics object to @@ -289,20 +278,19 @@ void InterfaceKinetics::_update_rates_C() * species for each phase. */ tp->getActivityConcentrations(DATA_PTR(m_actConc) + m_start[n]); - // - // Get regular concentrations too - // - tp->getConcentrations(DATA_PTR(m_conc) + m_start[n]); + + // Get regular concentrations too + tp->getConcentrations(DATA_PTR(m_conc) + m_start[n]); } m_ROP_ok = false; } -//============================================================================================================================ + void InterfaceKinetics::getActivityConcentrations(doublereal* const conc) { _update_rates_C(); copy(m_actConc.begin(), m_actConc.end(), conc); } -//============================================================================================================================ + void InterfaceKinetics::updateKc() { fill(m_rkcn.begin(), m_rkcn.end(), 0.0); @@ -312,9 +300,9 @@ void InterfaceKinetics::updateKc() * Get the vector of standard state electrochemical potentials for species in the Interfacial * kinetics object and store it in m_mu0[] and m_mu0_Kc[] */ - updateMu0(); - doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature()); - + updateMu0(); + doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature()); + // compute Delta mu^0 for all reversible reactions m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0_Kc), DATA_PTR(m_rkcn)); @@ -331,12 +319,10 @@ void InterfaceKinetics::updateKc() } } } -//============================================================================================================================ + void InterfaceKinetics::updateMu0() -{ - // +{ // First task is update the electrical potentials from the Phases - // _update_rates_phi(); updateExchangeCurrentQuantities(); @@ -348,21 +334,19 @@ void InterfaceKinetics::updateMu0() doublereal rt = GasConstant * thermo(0).temperature(); size_t np = nPhases(); for (size_t n = 0; n < np; n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); - nsp = thermo(n).nSpecies(); - for (size_t k = 0; k < nsp; k++) { - m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k); - m_mu0_Kc[ik] -= rt * thermo(n).logStandardConc(k); - ik++; - } + thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); + nsp = thermo(n).nSpecies(); + for (size_t k = 0; k < nsp; k++) { + m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k); + m_mu0_Kc[ik] -= rt * thermo(n).logStandardConc(k); + ik++; + } } } -//============================================================================================================================ + void InterfaceKinetics::checkPartialEquil() { - // // First task is update the electrical potentials from the Phases - // _update_rates_phi(); vector_fp dmu(nTotalSpecies(), 0.0); @@ -395,25 +379,25 @@ void InterfaceKinetics::checkPartialEquil() } } } -//============================================================================================================================ + void InterfaceKinetics::getFwdRatesOfProgress(doublereal* fwdROP) { updateROP(); std::copy(m_ropf.begin(), m_ropf.end(), fwdROP); } -//============================================================================================================================ + 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(); @@ -427,17 +411,17 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) kc[i] = exp(-kc[i]*rrt); } } -//=========================================================================================================== -/* - * values needed to convert from exchange current density to surface reaction rate. - * Calculate: - * - m_StandardConc[] - * - m_ProdStandConcReac[] - * - m_deltaG0[] - * - m_mu0[] - */ + void InterfaceKinetics::updateExchangeCurrentQuantities() { + /* + * Calculate: + * - m_StandardConc[] + * - m_ProdStandConcReac[] + * - m_deltaG0[] + * - m_mu0[] + */ + /* * First collect vectors of the standard Gibbs free energies of the * species and the standard concentrations @@ -457,38 +441,34 @@ 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(); m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net); } -//=========================================================================================================== + void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) { - // // Compute the electrical potential energy of each species - // size_t ik = 0; for (size_t n = 0; n < nPhases(); n++) { size_t nsp = thermo(n).nSpecies(); @@ -497,11 +477,10 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) ik++; } } - // + // Compute the change in electrical potential energy for each // reaction. This will only be non-zero if a potential // difference is present. - // m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot), DATA_PTR(deltaElectricEnergy_)); // Modify the reaction rates. Only modify those with a @@ -520,103 +499,83 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) #endif for (size_t i = 0; i < m_beta.size(); i++) { size_t irxn = m_ctrxn[i]; - // - // If we calculate the BV form directly, we don't add the voltage correction to the - // forward reaction rate constants. - // - if (m_ctrxn_BVform[i] == 0) { - eamod = m_beta[i] * deltaElectricEnergy_[irxn]; - if (eamod != 0.0) { + + // If we calculate the BV form directly, we don't add the voltage correction to the + // forward reaction rate constants. + if (m_ctrxn_BVform[i] == 0) { + eamod = m_beta[i] * deltaElectricEnergy_[irxn]; + if (eamod != 0.0) { #ifdef DEBUG_KIN_MODE - ea = GasConstant * m_E[irxn]; - if (eamod + ea < 0.0) { - writelog("Warning: act energy mod too large!\n"); - writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n"); - writelog(" Delta Ea = "+fp2str(eamod)+"\n"); - writelog(" Ea = "+fp2str(ea)+"\n"); - for (n = 0; n < np; n++) { - writelog("Phase "+int2str(n)+": phi = " - +fp2str(m_phi[n])+"\n"); - } - } + ea = GasConstant * m_E[irxn]; + if (eamod + ea < 0.0) { + writelog("Warning: act energy mod too large!\n"); + writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n"); + writelog(" Delta Ea = "+fp2str(eamod)+"\n"); + writelog(" Ea = "+fp2str(ea)+"\n"); + for (n = 0; n < np; n++) { + writelog("Phase "+int2str(n)+": phi = " + +fp2str(m_phi[n])+"\n"); + } + } #endif - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; - kf[irxn] *= exp(-eamod*rrt); - } - } + doublereal rt = GasConstant*thermo(0).temperature(); + doublereal rrt = 1.0/rt; + kf[irxn] *= exp(-eamod*rrt); + } + } } } -//================================================================================================================== -/* - * For a reaction rate constant that was given in units of Amps/m2 (exchange current - * density formulation with iECDFormulation == true), convert the rate to - * kmoles/m2/s. - * - * For a reaction rate constant that was given in units of kmol/m2/sec when the - * reaction type is a butler-volmer form, convert it to exchange current density - * form (amps/m2). - * - */ + void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* const kfwd) { updateExchangeCurrentQuantities(); doublereal rt = GasConstant * thermo(0).temperature(); doublereal rrt = 1.0/rt; - // + // Loop over all reactions which are defined to have a voltage transfer coefficient that // affects the activity energy for the reaction - // for (size_t i = 0; i < m_ctrxn.size(); i++) { - size_t irxn = m_ctrxn[i]; - // - // Determine whether the reaction rate constant is in an exchange current density formulation format. - // - int iECDFormulation = m_ctrxn_ecdf[i]; - if (iECDFormulation) { - // - // If the BV form is to be converted into the normal form then we go through this process. - // If it isn't to be converted, then we don't go through this process. - // - // We need to have the straight chemical reaction rate constant to come out of this calculation. - // - if (m_ctrxn_BVform[i] == 0) { - // - // Calculate the term and modify the forward reaction - // - double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt); - double tmp2 = m_ProdStanConcReac[irxn]; - tmp *= 1.0 / tmp2 / Faraday; - kfwd[irxn] *= tmp; - } - // - // If BVform is nonzero we don't need to do anything. - // - } else { - // - // kfwd[] is the chemical reaction rate constant - // - // If we are to calculate the BV form directly, then we will do the reverse. - // We will calculate the exchange current density formulation here and - // substitute it. - // - if (m_ctrxn_BVform[i] != 0) { - // - // Calculate the term and modify the forward reaction rate constant so that - // it's in the exchange current density formulation format - // - double tmp = exp(m_beta[i] * m_deltaG0[irxn] * rrt); - double tmp2 = m_ProdStanConcReac[irxn]; - tmp *= Faraday * tmp2; - kfwd[irxn] *= tmp; - } - } + size_t irxn = m_ctrxn[i]; + + // Determine whether the reaction rate constant is in an exchange current density formulation format. + int iECDFormulation = m_ctrxn_ecdf[i]; + if (iECDFormulation) { + // If the BV form is to be converted into the normal form then we go through this process. + // If it isn't to be converted, then we don't go through this process. + // + // We need to have the straight chemical reaction rate constant to come out of this calculation. + if (m_ctrxn_BVform[i] == 0) { + // + // Calculate the term and modify the forward reaction + // + double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt); + double tmp2 = m_ProdStanConcReac[irxn]; + tmp *= 1.0 / tmp2 / Faraday; + kfwd[irxn] *= tmp; + } + // If BVform is nonzero we don't need to do anything. + + } else { + // kfwd[] is the chemical reaction rate constant + // + // If we are to calculate the BV form directly, then we will do the reverse. + // We will calculate the exchange current density formulation here and + // substitute it. + if (m_ctrxn_BVform[i] != 0) { + + // Calculate the term and modify the forward reaction rate constant so that + // it's in the exchange current density formulation format + double tmp = exp(m_beta[i] * m_deltaG0[irxn] * rrt); + double tmp2 = m_ProdStanConcReac[irxn]; + tmp *= Faraday * tmp2; + kfwd[irxn] *= tmp; + } + } } } -//================================================================================================================== + void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) { - updateROP(); // copy rate coefficients into kfwd @@ -624,9 +583,8 @@ void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) // multiply by perturbation factor multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin()); - } -//================================================================================================================== + void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) { getFwdRateConstants(krev); @@ -639,95 +597,82 @@ void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversibl multiply_each(krev, krev + nReactions(), m_rkcn.begin()); } } -//============================================================================================================================ + void InterfaceKinetics::updateROP() { // evaluate rate constants and equilibrium constants at temperature and phi (electric potential) _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 the reaction rate coefficients, m_rfn, into m_ropf - // copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); - // + // Multiply by the perturbation factor - // multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin()); // // 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 // 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. - // + // 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; - + 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; - } - /* + getDeltaGibbs(0); - double exp1 = nu * nStoich * beta / rtdf - double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (GasConstant * temp); - double val = io * (exp(exp1) - exp(exp2)); + if (nStoichElectrons != 0.0) { + OCV = m_deltaG[jrxn]/Faraday/ nStoichElectrons; + } - 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]; - */ + /* + + 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 (size_t j = 0; j != m_ii; ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; } @@ -795,32 +740,30 @@ void InterfaceKinetics::updateROP() m_ROP_ok = true; } -//================================================================================================================== + void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) { /* - * Get the chemical potentials of the species in the all of the phases used in the + * Get the chemical potentials of the species in the all of the phases used in the * kinetics mechanism */ for (size_t n = 0; n < nPhases(); n++) { m_thermo[n]->getChemPotentials(DATA_PTR(m_mu) + m_start[n]); } - // - // Use the stoichiometric manager to find deltaG for each - // reaction. - // + + // 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]; - } + for (size_t j = 0; j < m_ii; ++j) { + deltaG[j] = m_deltaG[j]; + } } } -//================================================================================================================== + void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) { /* - * Get the chemical potentials of the species + * Get the chemical potentials of the species */ size_t np = nPhases(); for (size_t n = 0; n < np; n++) { @@ -832,7 +775,7 @@ void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaM); } -//================================================================================================================== + void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) { /* @@ -847,7 +790,7 @@ void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); } -//=========================================================================================================== + void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS) { /* @@ -863,7 +806,7 @@ void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } -//=========================================================================================================== + void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS) { /* @@ -881,7 +824,7 @@ void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), deltaGSS); } -//=========================================================================================================== + void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) { /* @@ -903,7 +846,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); } -//=========================================================================================================== + void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS) { /* @@ -924,7 +867,7 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } -//============================================================================================================================ + void InterfaceKinetics::addReaction(ReactionData& r) { int reactionType = r.reactionType; @@ -932,19 +875,18 @@ void InterfaceKinetics::addReaction(ReactionData& r) reactionType_.push_back(reactionType); if ((reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN ) || - (reactionType == BUTLERVOLMER_RXN ) || - (reactionType == SURFACEAFFINITY_RXN) || - (reactionType == GLOBAL_RXN)) { - // - // Add global reactions - // - addGlobalReaction(r); + (reactionType == BUTLERVOLMER_RXN ) || + (reactionType == SURFACEAFFINITY_RXN) || + (reactionType == GLOBAL_RXN)) { + + // Add global reactions + addGlobalReaction(r); } else { - /* - * Install the rate coefficient for the current reaction - * in the appropriate data structure. - */ - addElementaryReaction(r); + /* + * Install the rate coefficient for the current reaction + * in the appropriate data structure. + */ + addElementaryReaction(r); } /* * Add the reactants and products for m_ropnet;the current reaction @@ -963,7 +905,7 @@ void InterfaceKinetics::addReaction(ReactionData& r) * increase the size of m_perturb by one as well. */ incrementRxnCount(); - // + m_rxneqn.push_back(r.equation); m_rxnPhaseIsReactant.push_back(std::vector(nPhases(), false)); @@ -983,27 +925,22 @@ void InterfaceKinetics::addReaction(ReactionData& r) m_rxnPhaseIsProduct[i][p] = true; } } -//============================================================================================================================ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) { - // // install rate coefficient calculator - // vector_fp& rp = rdata.rateCoeffParameters; size_t ncov = rdata.cov.size(); - // + // Turn on the global flag indicating surface coverage dependence - // if (ncov > 3) { m_has_coverage_dependence = true; } for (size_t m = 0; m < ncov; m++) { rp.push_back(rdata.cov[m]); } - // + // Find out the reaction type - // int reactionType = rdata.reactionType; reactionTypes_.push_back(reactionType); @@ -1035,33 +972,31 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) m_has_electrochem_rxns = true; m_beta.push_back(rdata.beta); m_ctrxn.push_back(m_ii); - m_ctrxn_BVform.push_back(0); + m_ctrxn_BVform.push_back(0); if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) { m_has_exchange_current_density_formulation = true; m_ctrxn_ecdf.push_back(1); } else { m_ctrxn_ecdf.push_back(0); } - m_ctrxn_ROPOrdersList_.push_back(0); - m_ctrxn_FwdOrdersList_.push_back(0); - if (rdata.filmResistivity > 0.0) { - throw CanteraError("InterfaceKinetics::addElementaryReaction()", - "film resistivity set for elementary reaction"); - } - m_ctrxn_resistivity_.push_back(rdata.filmResistivity); + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); + if (rdata.filmResistivity > 0.0) { + throw CanteraError("InterfaceKinetics::addElementaryReaction()", + "film resistivity set for elementary reaction"); + } + m_ctrxn_resistivity_.push_back(rdata.filmResistivity); } // add constant term to rate coeff value vector m_rfn.push_back(rdata.rateCoeffParameters[0]); registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc); } -//============================================================================================================================ + void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) { - // // Install rate coeff calculator // This is done no matter what the type of reaction it is - // vector_fp& rp = rdata.rateCoeffParameters; size_t ncov = rdata.cov.size(); if (ncov > 3) { @@ -1071,9 +1006,7 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) rp.push_back(rdata.cov[m]); } - // // Find out the reaction type - // int reactionType = rdata.reactionType; reactionTypes_.push_back(reactionType); @@ -1097,29 +1030,25 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) * Change the reaction rate coefficient type back to its original value */ rdata.rateCoeffType = reactionRateCoeffType_orig; - // + // Store activation energy - // m_E.push_back(rdata.rateCoeffParameters[2]); - // // Add the reaction into the list of electrochemical extras - // if (rdata.beta > 0.0 || 1) { m_has_electrochem_rxns = true; m_beta.push_back(rdata.beta); - // Push back the id of the reaction + // Push back the id of the reaction m_ctrxn.push_back(m_ii); - // + // Specify alternative forms of the electrochemical reaction - // if (rdata.reactionType == BUTLERVOLMER_RXN) { - m_ctrxn_BVform.push_back(1); + m_ctrxn_BVform.push_back(1); } else if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) { - m_ctrxn_BVform.push_back(2); + m_ctrxn_BVform.push_back(2); } else { - // set the default to be the normal forward / reverse calculation method - m_ctrxn_BVform.push_back(0); + // 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; @@ -1127,38 +1056,35 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) } else { m_ctrxn_ecdf.push_back(0); } - // + // Store the film resistivity - // m_ctrxn_resistivity_.push_back(rdata.filmResistivity); - - if (rdata.forwardFullOrder_.size() > 0) { - RxnOrders* ro = new RxnOrders(); - ro->fill(rdata.forwardFullOrder_); - m_ctrxn_ROPOrdersList_.push_back(ro); - m_ctrxn_FwdOrdersList_.push_back(0); - // - // - // Fill in the Fwd Orders dependence here for B-V reactions - // - if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) { - std::vector fwdFullorders(m_kk, 0.0); - determineFwdOrdersBV(rdata, fwdFullorders); - RxnOrders* ro = new RxnOrders(); - ro->fill(fwdFullorders); - m_ctrxn_FwdOrdersList_[m_ii] = ro; - } - } else { - m_ctrxn_ROPOrdersList_.push_back(0); - m_ctrxn_FwdOrdersList_.push_back(0); - } + + if (rdata.forwardFullOrder_.size() > 0) { + RxnOrders* ro = new RxnOrders(); + ro->fill(rdata.forwardFullOrder_); + m_ctrxn_ROPOrdersList_.push_back(ro); + m_ctrxn_FwdOrdersList_.push_back(0); + + // Fill in the Fwd Orders dependence here for B-V reactions + if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) { + std::vector fwdFullorders(m_kk, 0.0); + determineFwdOrdersBV(rdata, fwdFullorders); + RxnOrders* ro = new RxnOrders(); + ro->fill(fwdFullorders); + m_ctrxn_FwdOrdersList_[m_ii] = ro; + } + } else { + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); + } } // add constant term to rate coeff value vector m_rfn.push_back(rdata.rateCoeffParameters[0]); registerReaction(m_ii, ELEMENTARY_RXN, iloc); } -//================================================================================================================== + void InterfaceKinetics::setIOFlag(int ioFlag) { m_ioFlag = ioFlag; @@ -1166,10 +1092,9 @@ void InterfaceKinetics::setIOFlag(int ioFlag) m_integrator->setIOFlag(ioFlag); } } -//================================================================================================================== + void InterfaceKinetics::installReagents(const ReactionData& r) { - size_t n, ns, m; doublereal nsFlt; /* @@ -1260,14 +1185,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); } -//================================================================================================================== + void InterfaceKinetics::init() { m_kk = 0; @@ -1285,7 +1210,7 @@ void InterfaceKinetics::init() m_pot.resize(m_kk, 0.0); m_phi.resize(nPhases(), 0.0); } -//================================================================================================================== + void InterfaceKinetics::finalize() { Kinetics::finalize(); @@ -1294,9 +1219,8 @@ void InterfaceKinetics::finalize() size_t ks = reactionPhaseIndex(); if (ks == npos) throw CanteraError("InterfaceKinetics::finalize", "no surface phase is present."); - // + // Check to see that the interface routine has a dimension of 2 - // m_surf = (SurfPhase*)&thermo(ks); if (m_surf->nDim() != 2) { throw CanteraError("InterfaceKinetics::finalize", @@ -1322,9 +1246,8 @@ void InterfaceKinetics::finalize() m_ropnet.resize(1, 0.0); m_rkcn.resize(1, 0.0); } - // + // Malloc and calculate all of the quantities that go into the extra description of reactions - // rmcVector.resize(m_ii, 0); for (size_t i = 0; i < m_ii; i++) { rmcVector[i] = new RxnMolChange(this, i); @@ -1332,7 +1255,7 @@ void InterfaceKinetics::finalize() m_finalized = true; } -//================================================================================================================== + doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const { for (size_t i = 0; i < m_ctrxn.size(); i++) { @@ -1342,12 +1265,12 @@ doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const } return 0.0; } -//================================================================================================================== + bool InterfaceKinetics::ready() const { return m_finalized; } -//================================================================================================================== + void InterfaceKinetics::advanceCoverages(doublereal tstep) { if (m_integrator == 0) { @@ -1360,7 +1283,7 @@ void InterfaceKinetics::advanceCoverages(doublereal tstep) delete m_integrator; m_integrator = 0; } -//================================================================================================================== + void InterfaceKinetics::solvePseudoSteadyStateProblem( int ifuncOverride, doublereal timeScaleOverride) { @@ -1377,7 +1300,7 @@ void InterfaceKinetics::solvePseudoSteadyStateProblem( */ m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride); } -//================================================================================================================== + void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) { if (iphase >= m_thermo.size()) { @@ -1397,9 +1320,8 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) } m_phaseIsStable[iphase] = false; } - } -//================================================================================================================== + int InterfaceKinetics::phaseExistence(const size_t iphase) const { if (iphase >= m_thermo.size()) { @@ -1407,7 +1329,7 @@ int InterfaceKinetics::phaseExistence(const size_t iphase) const } return m_phaseExists[iphase]; } -//================================================================================================================== + int InterfaceKinetics::phaseStability(const size_t iphase) const { if (iphase >= m_thermo.size()) { @@ -1415,22 +1337,22 @@ int InterfaceKinetics::phaseStability(const size_t iphase) const } return m_phaseIsStable[iphase]; } -//================================================================================================================== + doublereal InterfaceKinetics::reactantStoichCoeff(size_t kSpecKin, size_t irxn) const { return getValue(m_rrxn[kSpecKin], irxn, 0.0); -} -//================================================================================================================== -doublereal InterfaceKinetics::productStoichCoeff(size_t kSpecKin, size_t irxn) const +} + +doublereal InterfaceKinetics::productStoichCoeff(size_t kSpecKin, size_t irxn) const { return getValue(m_prxn[kSpecKin], irxn, 0.0); } -//================================================================================================================== -int InterfaceKinetics::reactionType(size_t irxn) const + +int InterfaceKinetics::reactionType(size_t irxn) const { return reactionType_[irxn]; } -//================================================================================================================== + void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable) { if (iphase >= m_thermo.size()) { @@ -1442,68 +1364,52 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl m_phaseIsStable[iphase] = false; } } -//================================================================================================================== -// 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) + +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 InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders) { - // // Start out with the full ROP orders vector. // This vector will have the BV exchange current density orders in it. - // fwdFullorders = rdata.forwardFullOrder_; - // + // forward and reverse beta values - // double betaf = rdata.beta; double betar = 1.0 - betaf; - // + // Loop over the reactants doing away with the BV terms. // This should leave the reactant terms only, even if they are non-mass action. - // for (size_t j = 0; j < rdata.reactants.size(); j++) { - size_t kkin = rdata.reactants[j]; - double oo = rdata.rstoich[j]; - fwdFullorders[kkin] += betaf * oo; - // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet) - if (abs(fwdFullorders[kkin]) < 0.00001) { - fwdFullorders[kkin] = 0.0; - } + size_t kkin = rdata.reactants[j]; + double oo = rdata.rstoich[j]; + fwdFullorders[kkin] += betaf * oo; + // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet) + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } } - // + // Loop over the products doing away with the BV terms. // This should leave the reactant terms only, even if they are non-mass action. - // for (size_t j = 0; j < rdata.products.size(); j++) { - size_t kkin = rdata.products[j]; - double oo = rdata.pstoich[j]; - fwdFullorders[kkin] -= betaf * oo; - if (abs(fwdFullorders[kkin]) < 0.00001) { - fwdFullorders[kkin] = 0.0; - } + size_t kkin = rdata.products[j]; + double oo = rdata.pstoich[j]; + fwdFullorders[kkin] -= betaf * oo; + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } } } -//================================================================================================================== + void EdgeKinetics::finalize() { - // // Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below. - // Therefore, we have to malloc room in arrays that would normally be + // Therefore, we have to malloc room in arrays that would normally be // handled by the InterfaceKinetics::finalize() call. - // Kinetics::finalize(); size_t safe_reaction_size = std::max(m_ii, 1); @@ -1511,9 +1417,8 @@ void EdgeKinetics::finalize() size_t ks = reactionPhaseIndex(); if (ks == npos) throw CanteraError("EdgeKinetics::finalize", "no surface phase is present."); - // + // Check to see edge phase has a dimension of 1 - // m_surf = (SurfPhase*)&thermo(ks); if (m_surf->nDim() != 1) { throw CanteraError("EdgeKinetics::finalize", @@ -1539,9 +1444,8 @@ void EdgeKinetics::finalize() m_ropnet.resize(1, 0.0); m_rkcn.resize(1, 0.0); } - // + // Malloc and calculate all of the quantities that go into the extra description of reactions - // rmcVector.resize(m_ii, 0); for (size_t i = 0; i < m_ii; i++) { rmcVector[i] = new RxnMolChange(this, i); @@ -1549,21 +1453,21 @@ void EdgeKinetics::finalize() m_finalized = true; } -//================================================================================================================== + RxnOrders::RxnOrders() { } -//================================================================================================================== + RxnOrders::~RxnOrders() { } -//================================================================================================================== + RxnOrders::RxnOrders(const RxnOrders& right) : kinSpeciesIDs_(right.kinSpeciesIDs_), - kinSpeciesOrders_(right.kinSpeciesOrders_) + kinSpeciesOrders_(right.kinSpeciesOrders_) { } -//================================================================================================================== + RxnOrders& RxnOrders::operator=(const RxnOrders& right) { if (this == &right) { @@ -1573,20 +1477,20 @@ RxnOrders& RxnOrders::operator=(const RxnOrders& right) kinSpeciesOrders_ = right.kinSpeciesOrders_; return *this; } -//================================================================================================================== + int RxnOrders::fill(const std::vector& fullForwardOrders) { int nzeroes = 0; kinSpeciesIDs_.clear(); kinSpeciesOrders_.clear(); - for (size_t k = 0; k < fullForwardOrders.size(); ++k) { - if (fullForwardOrders[k] != 0.0) { - kinSpeciesIDs_.push_back(k); - kinSpeciesOrders_.push_back(fullForwardOrders[k]); + for (size_t k = 0; k < fullForwardOrders.size(); ++k) { + if (fullForwardOrders[k] != 0.0) { + kinSpeciesIDs_.push_back(k); + kinSpeciesOrders_.push_back(fullForwardOrders[k]); ++nzeroes; - } + } } return nzeroes; } -//================================================================================================================== + }