diff --git a/Cantera/src/kinetics/AqueousKinetics.cpp b/Cantera/src/kinetics/AqueousKinetics.cpp new file mode 100644 index 000000000..0559dec2b --- /dev/null +++ b/Cantera/src/kinetics/AqueousKinetics.cpp @@ -0,0 +1,550 @@ +/** + * @file AqueousKinetics.cpp + * + * Homogeneous kinetics in an aqueous phase, either condensed + * or dilute in salts + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Date$ + * $Revision$ + */ + + + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "AqueousKinetics.h" +#include "ReactionData.h" +#include "RateCoeffMgr.h" + +#include +using namespace std; + +namespace Cantera { + + /** + * Construct an empty reaction mechanism. + */ + AqueousKinetics::AqueousKinetics(thermo_t* thermo) : + Kinetics(), + m_kk(0), + m_nfall(0), + m_nirrev(0), + m_nrev(0), + m_finalized(false) + { + if (thermo != 0) addPhase(*thermo); + m_kdata = new AqueousKineticsData; + m_kdata->m_temp = 0.0; + m_rxnstoich = new ReactionStoichMgr; + } + + AqueousKinetics::~AqueousKinetics() { + delete m_kdata; + delete m_rxnstoich; + } + + /** + * Update temperature-dependent portions of reaction rates and + * falloff functions. + */ + void AqueousKinetics:: + update_T() {} + + void AqueousKinetics:: + update_C() {} + + void AqueousKinetics::_update_rates_T() { + doublereal T = thermo().temperature(); + // m_kdata->m_logStandConc = log(thermo().standardConcentration()); + doublereal logT = log(T); + m_rates.update(T, logT, &m_kdata->m_rfn[0]); + + + m_kdata->m_temp = T; + updateKc(); + m_kdata->m_ROP_ok = false; + + }; + + + /** + * Update properties that depend on concentrations. Currently only + * the enhanced collision partner concentrations are updated here. + */ + void AqueousKinetics:: + _update_rates_C() { + thermo().getActivityConcentrations(&m_conc[0]); + + m_kdata->m_ROP_ok = false; + } + + /** + * Update the equilibrium constants in molar units. + */ + void AqueousKinetics::updateKc() { + int i, irxn; + vector_fp& m_rkc = m_kdata->m_rkcn; + doublereal rt = GasConstant* m_kdata->m_temp; + + thermo().getStandardChemPotentials(&m_grt[0]); + fill(m_rkc.begin(), m_rkc.end(), 0.0); + int nsp = thermo().nSpecies(); + for (int k = 0; k < nsp; k++) { + doublereal logStandConc_k = thermo().logStandardConc(k); + m_grt[k] -= rt * logStandConc_k; + } + + // compute Delta G^0 for all reversible reactions + m_rxnstoich->getRevReactionDelta(m_ii, &m_grt[0], &m_rkc[0]); + + //doublereal logStandConc = m_kdata->m_logStandConc; + doublereal rrt = 1.0/(GasConstant * thermo().temperature()); + for (i = 0; i < m_nrev; i++) { + irxn = m_revindex[i]; + m_rkc[irxn] = exp(m_rkc[irxn]*rrt); + } + + for(i = 0; i != m_nirrev; ++i) { + m_rkc[ m_irrev[i] ] = 0.0; + } + } + + /** + * Get the equilibrium constants of all reactions, whether + * reversible or not. + */ + void AqueousKinetics::getEquilibriumConstants(doublereal* kc) { + int i; + _update_rates_T(); + vector_fp& rkc = m_kdata->m_rkcn; + + thermo().getStandardChemPotentials(&m_grt[0]); + fill(rkc.begin(), rkc.end(), 0.0); + doublereal rt = GasConstant * m_kdata->m_temp; + int nsp = thermo().nSpecies(); + for (int k = 0; k < nsp; k++) { + doublereal logStandConc_k = thermo().logStandardConc(k); + m_grt[k] -= rt * logStandConc_k; + } + + // compute Delta G^0 for all reactions + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], &rkc[0]); + + doublereal rrt = 1.0/(GasConstant * thermo().temperature()); + for (i = 0; i < m_ii; i++) { + kc[i] = exp(-rkc[i]*rrt); + } + + // force an update of T-dependent properties, so that m_rkcn will + // be updated before it is used next. + m_kdata->m_temp = 0.0; + } + + /** + * + * getDeltaGibbs(): + * + * Return the vector of values for the reaction gibbs free energy + * change + * These values depend upon the concentration + * of the ideal gas. + * + * units = J kmol-1 + */ + void AqueousKinetics::getDeltaGibbs(doublereal* deltaG) { + /* + * Get the chemical potentials of the species in the + * ideal gas solution. + */ + thermo().getChemPotentials(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaG); + } + + /** + * + * getDeltaEnthalpy(): + * + * Return the vector of values for the reactions change in + * enthalpy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH) { + /* + * Get the partial molar enthalpy of all species in the + * ideal gas. + */ + thermo().getPartialMolarEnthalpies(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaH); + } + + /* + * + * getDeltaEntropy(): + * + * Return the vector of values for the reactions change in + * entropy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + void AqueousKinetics::getDeltaEntropy( doublereal* deltaS) { + /* + * Get the partial molar entropy of all species in the + * solid solution. + */ + thermo().getPartialMolarEntropies(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaS for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaS); + } + + /** + * + * getDeltaSSGibbs(): + * + * Return the vector of values for the reaction + * standard state gibbs free energy change. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG) { + /* + * Get the standard state chemical potentials of the species. + * This is the array of chemical potentials at unit activity + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getStandardChemPotentials(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaG); + } + + /** + * + * getDeltaSSEnthalpy(): + * + * Return the vector of values for the change in the + * standard state enthalpies of reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH) { + /* + * Get the standard state enthalpies of the species. + * This is the array of chemical potentials at unit activity + * We define these here as the enthalpies of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getEnthalpy_RT(&m_grt[0]); + doublereal RT = thermo().temperature() * GasConstant; + for (int k = 0; k < m_kk; k++) { + m_grt[k] *= RT; + } + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaH); + } + + /* + * + * getDeltaSSEntropy(): + * + * Return the vector of values for the change in the + * standard state entropies for each reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS) { + /* + * Get the standard state entropy of the species. + * We define these here as the entropies of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getEntropy_R(&m_grt[0]); + doublereal R = GasConstant; + for (int k = 0; k < m_kk; k++) { + m_grt[k] *= R; + } + /* + * Use the stoichiometric manager to find deltaS for each + * reaction. + */ + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaS); + } + + + + void AqueousKinetics::updateROP() { + + _update_rates_T(); + _update_rates_C(); + + if (m_kdata->m_ROP_ok) return; + + const vector_fp& rf = m_kdata->m_rfn; + const vector_fp& m_rkc = m_kdata->m_rkcn; + array_fp& ropf = m_kdata->m_ropf; + array_fp& ropr = m_kdata->m_ropr; + array_fp& ropnet = m_kdata->m_ropnet; + + // copy rate coefficients into ropf + copy(rf.begin(), rf.end(), ropf.begin()); + + + // multiply by perturbation factor + multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); + + // copy the forward rates to the reverse rates + copy(ropf.begin(), ropf.end(), 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(ropr.begin(), ropr.end(), m_rkc.begin()); + + // multiply ropf by concentration products + m_rxnstoich->multiplyReactants(&m_conc[0], &ropf[0]); + //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); + + // for reversible reactions, multiply ropr by concentration + // products + m_rxnstoich->multiplyRevProducts(&m_conc[0], &ropr[0]); + //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); + + for (int j = 0; j != m_ii; ++j) { + ropnet[j] = ropf[j] - ropr[j]; + } + + m_kdata->m_ROP_ok = true; + } + + /** + * + * getFwdRateConstants(): + * + * Update the rate of progress for the reactions. + * This key routine makes sure that the rate of progress vectors + * located in the solid kinetics data class are up to date. + */ + void AqueousKinetics:: + getFwdRateConstants(doublereal *kfwd) { + _update_rates_T(); + _update_rates_C(); + + // copy rate coefficients into ropf + const vector_fp& rf = m_kdata->m_rfn; + array_fp& ropf = m_kdata->m_ropf; + copy(rf.begin(), rf.end(), ropf.begin()); + + + + // multiply by perturbation factor + multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); + + for (int i = 0; i < m_ii; i++) { + kfwd[i] = ropf[i]; + } + } + + /** + * + * getRevRateConstants(): + * + * Return a vector of the reverse reaction rate constants + * + * Length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. + */ + void AqueousKinetics:: + getRevRateConstants(doublereal *krev, bool doIrreversible) { + /* + * go get the forward rate constants. -> note, we don't + * really care about speed or redundancy in these + * informational routines. + */ + getFwdRateConstants(krev); + + if (doIrreversible) { + doublereal *tmpKc = &m_kdata->m_ropnet[0]; + getEquilibriumConstants(tmpKc); + for (int i = 0; i < m_ii; i++) { + krev[i] /= tmpKc[i]; + } + } else { + /* + * m_rkc[] is zero for irreversibly reactions + */ + const vector_fp& m_rkc = m_kdata->m_rkcn; + for (int i = 0; i < m_ii; i++) { + krev[i] *= m_rkc[i]; + } + } + } + + void AqueousKinetics::addReaction(const ReactionData& r) { + + if (r.reactionType == ELEMENTARY_RXN) addElementaryReaction(r); + + // operations common to all reaction types + installReagents( r ); + installGroups(reactionNumber(), r.rgroups, r.pgroups); + incrementRxnCount(); + m_rxneqn.push_back(r.equation); + } + + + + + void AqueousKinetics::addElementaryReaction(const ReactionData& r) { + int iloc; + + // install rate coeff calculator + iloc = m_rates.install( reactionNumber(), + r.rateCoeffType, r.rateCoeffParameters.size(), + DATA_PTR(r.rateCoeffParameters) ); + + // add constant term to rate coeff value vector + m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); + + // forward rxn order equals number of reactants + m_fwdOrder.push_back(r.reactants.size()); + registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc); + } + + + + + void AqueousKinetics::installReagents(const ReactionData& r) { + + m_kdata->m_ropf.push_back(0.0); // extend by one for new rxn + m_kdata->m_ropr.push_back(0.0); + m_kdata->m_ropnet.push_back(0.0); + int n, ns, m; + doublereal nsFlt; + doublereal reactantGlobalOrder = 0.0; + doublereal productGlobalOrder = 0.0; + int rnum = reactionNumber(); + + vector_int rk; + int nr = r.reactants.size(); + for (n = 0; n < nr; n++) { + nsFlt = r.rstoich[n]; + reactantGlobalOrder += nsFlt; + ns = (int) nsFlt; + if ((doublereal) ns != nsFlt) { + if (ns < 1) { + ns = 1; + } + } + if (r.rstoich[n] != 0.0) + m_rrxn[r.reactants[n]][rnum] += r.rstoich[n]; + for (m = 0; m < ns; m++) { + rk.push_back(r.reactants[n]); + } + } + m_reactants.push_back(rk); + + vector_int pk; + int np = r.products.size(); + for (n = 0; n < np; n++) { + nsFlt = r.pstoich[n]; + productGlobalOrder += nsFlt; + ns = (int) nsFlt; + if ((double) ns != nsFlt) { + if (ns < 1) { + ns = 1; + } + } + if (r.pstoich[n] != 0.0) + m_prxn[r.products[n]][rnum] += r.pstoich[n]; + for (m = 0; m < ns; m++) { + pk.push_back(r.products[n]); + } + } + m_products.push_back(pk); + + m_kdata->m_rkcn.push_back(0.0); + + m_rxnstoich->add(reactionNumber(), r); + + if (r.reversible) { + m_dn.push_back(productGlobalOrder - reactantGlobalOrder); + m_revindex.push_back(reactionNumber()); + m_nrev++; + } + else { + m_dn.push_back(productGlobalOrder - reactantGlobalOrder); + m_irrev.push_back( reactionNumber() ); + m_nirrev++; + } + } + + + void AqueousKinetics::installGroups(int irxn, + const vector& r, + const vector& p) { + if (!r.empty()) { + writelog("installing groups for reaction "+int2str(reactionNumber())); + m_rgroups[reactionNumber()] = r; + m_pgroups[reactionNumber()] = p; + } + } + + + void AqueousKinetics::init() { + m_kk = thermo().nSpecies(); + m_rrxn.resize(m_kk); + m_prxn.resize(m_kk); + m_conc.resize(m_kk); + m_grt.resize(m_kk); + m_kdata->m_logp_ref = log(thermo().refPressure()) - log(GasConstant); + } + + void AqueousKinetics::finalize() { + if (!m_finalized) { + m_finalized = true; + } + } + + bool AqueousKinetics::ready() const { + return (m_finalized); + } + +} diff --git a/Cantera/src/kinetics/AqueousKinetics.h b/Cantera/src/kinetics/AqueousKinetics.h new file mode 100644 index 000000000..a23bc0626 --- /dev/null +++ b/Cantera/src/kinetics/AqueousKinetics.h @@ -0,0 +1,405 @@ +/** + * @file AqueousKinetics.h + * + * @ingroup chemkinetics + * + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_AQUEOUSKINETICS_H +#define CT_AQUEOUSKINETICS_H + +#include +#include + +#include "mix_defs.h" +#include "Kinetics.h" + +#include "utilities.h" + +#include "ReactionStoichMgr.h" +#include "ThirdBodyMgr.h" +#include "FalloffMgr.h" +#include "RateCoeffMgr.h" + +#include +#include + +void get_wdot(const doublereal* rop, doublereal* wdot); + +namespace Cantera { + + // forward references + + + class ReactionData; + class AqueousKineticsData; + class Thermo; + + /** + * Holds mechanism-specific data. + */ + class AqueousKineticsData { + public: + AqueousKineticsData() : + m_logp_ref(0.0), + m_logc_ref(0.0), + m_ROP_ok(false), + m_temp(0.0) + {} + virtual ~AqueousKineticsData(){} + + doublereal m_logp_ref, m_logc_ref; + array_fp m_ropf; + array_fp m_ropr, m_ropnet; + array_fp m_rfn_low, m_rfn_high; + bool m_ROP_ok; + + doublereal m_temp; + array_fp m_rfn; + + array_fp m_rkcn; + }; + + + /** + * Kinetics manager for elementary aqueous-phase chemistry. This + * kinetics manager implements standard mass-action reaction rate + * expressions for liquids + * + * + * Concentration + * + * @ingroup kinetics + */ + class AqueousKinetics : public Kinetics { + + public: + + /** + * @name Constructors and General Information + */ + //@{ + /// Constructor. + AqueousKinetics(thermo_t* thermo = 0); + + /// Destructor. + virtual ~AqueousKinetics(); + + virtual int ID() const { return cAqueousKinetics; } + virtual int type() const { return cAqueousKinetics; } + + virtual doublereal reactantStoichCoeff(int k, int i) const { + return m_rrxn[k][i]; + } + + virtual doublereal productStoichCoeff(int k, int i) const { + return m_prxn[k][i]; + } + + //@} + /** + * @name Reaction Rates Of Progress + */ + //@{ + /** + * Forward rates of progress. + * Return the forward rates of progress in array fwdROP, which + * must be dimensioned at least as large as the total number + * of reactions. + */ + virtual void getFwdRatesOfProgress(doublereal* fwdROP) { + updateROP(); + std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP); + } + + /** + * Reverse rates of progress. + * Return the reverse rates of progress in array revROP, which + * must be dimensioned at least as large as the total number + * of reactions. + */ + virtual void getRevRatesOfProgress(doublereal* revROP) { + updateROP(); + std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP); + } + + /** + * Net rates of progress. Return the net (forward - reverse) + * rates of progress in array netROP, which must be + * dimensioned at least as large as the total number of + * reactions. + */ + virtual void getNetRatesOfProgress(doublereal* netROP) { + updateROP(); + std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); + } + + + /** + * Equilibrium constants. Return the equilibrium constants of + * the reactions in concentration units in array kc, which + * must be dimensioned at least as large as the total number + * of reactions. + */ + virtual void getEquilibriumConstants(doublereal* kc); + + /** + * Return the array of values for the reaction gibbs free energy + * change. + * These values depend on the species concentrations. + * + * units = J kmol-1 + */ + virtual void getDeltaGibbs( doublereal* deltaG); + + /** + * Return the array of values for the reaction enthalpy change. + * These values depend upon the species concentrations. + * + * units = J kmol-1 + */ + virtual void getDeltaEnthalpy( doublereal* deltaH); + + /** + * Return the array of values for the reactions change in + * entropy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + virtual void getDeltaEntropy(doublereal* deltaS); + + /** + * Return the array of values for the reaction + * standard state Gibbs free energy change. + * These values do not depend on the species + * concentrations. + * + * units = J kmol-1 + */ + virtual void getDeltaSSGibbs(doublereal* deltaG); + + /** + * Return the array of values for the change in the + * standard state enthalpies of reaction. + * These values do not depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaSSEnthalpy(doublereal* deltaH); + + /** + * Return the array of values for the change in the + * standard state entropies for each reaction. + * These values do not depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + virtual void getDeltaSSEntropy(doublereal* deltaS); + + //@} + /** + * @name Species Production Rates + */ + //@{ + + //! Return the species net production rates + /*! + * Species net production rates [kmol/m^3/s]. Return the species + * net production rates (creation - destruction) in array + * wdot, which must be dimensioned at least as large as the + * total number of species. + * + * @param net Array of species production rates. + * units kmol m-3 s-1 + */ + virtual void getNetProductionRates(doublereal* net) { + updateROP(); + //#ifdef HWMECH + //get_wdot(&m_kdata->m_ropnet[0], net); + //#else + m_rxnstoich->getNetProductionRates(m_kk, + &m_kdata->m_ropnet[0], net); + //#endif + } + + /** + * Species creation rates [kmol/m^3]. Return the species + * creation rates in array cdot, which must be + * dimensioned at least as large as the total number of + * species. + * + */ + virtual void getCreationRates(doublereal* cdot) { + updateROP(); + m_rxnstoich->getCreationRates(m_kk, &m_kdata->m_ropf[0], + &m_kdata->m_ropr[0], cdot); + } + + /** + * Species destruction rates [kmol/m^3]. Return the species + * destruction rates in array ddot, which must be + * dimensioned at least as large as the total number of + * species. + * + */ + virtual void getDestructionRates(doublereal* ddot) { + updateROP(); + m_rxnstoich->getDestructionRates(m_kk, &m_kdata->m_ropf[0], + &m_kdata->m_ropr[0], ddot); + + } + + //@} + /** + * @name Reaction Mechanism Informational Query Routines + */ + //@{ + + /** + * Flag specifying the type of reaction. The legal values and + * their meaning are specific to the particular kinetics + * manager. + */ + virtual int reactionType(int i) const { + return m_index[i].first; + } + + virtual std::string reactionString(int i) const { + return m_rxneqn[i]; + } + + /** + * True if reaction i has been declared to be reversible. If + * isReversible(i) is false, then the reverse rate of progress + * for reaction i is always zero. + */ + virtual bool isReversible(int i) { + if (std::find(m_revindex.begin(), m_revindex.end(), i) + < m_revindex.end()) return true; + else return false; + } + + /** + * Return the forward rate constants + * + * length is the number of reactions. units depends + * on many issues. + */ + virtual void getFwdRateConstants(doublereal *kfwd); + + /** + * Return the reverse rate constants. + * + * length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. + */ + virtual void getRevRateConstants(doublereal *krev, + bool doIrreversible = false); + + //@} + /** + * @name Reaction Mechanism Setup Routines + */ + //@{ + + virtual void init(); + + /// Add a reaction to the mechanism. + virtual void addReaction(const ReactionData& r); + + virtual void finalize(); + virtual bool ready() const; + + virtual void update_T(); + virtual void update_C(); + + void updateROP(); + + + const std::vector& reactantGroups(int i) + { return m_rgroups[i]; } + const std::vector& productGroups(int i) + { return m_pgroups[i]; } + + + void _update_rates_T(); + void _update_rates_C(); + + //@} + + protected: + + int m_kk, m_nfall; + + Rate1 m_rates; + + mutable std::map > m_index; + + std::vector m_irrev; + + ReactionStoichMgr* m_rxnstoich; + + std::vector m_fwdOrder; + + int m_nirrev; + int m_nrev; + + std::map > m_rgroups; + std::map > m_pgroups; + + std::vector m_rxntype; + + mutable std::vector > m_rrxn; + mutable std::vector > m_prxn; + + /** + * Difference between the input global reactants order + * and the input global products order. Changed to a double + * to account for the fact that we can have real-valued + * stoichiometries. + */ + array_fp m_dn; + array_int m_revindex; + + std::vector m_rxneqn; + + AqueousKineticsData* m_kdata; + + array_fp m_conc; + array_fp m_grt; + + + private: + + int reactionNumber(){ return m_ii;} + std::vector > m_stoich; + + void addElementaryReaction(const ReactionData& r); + + + void installReagents(const ReactionData& r); + + void installGroups(int irxn, const std::vector& r, + const std::vector& p); + void updateKc(); + + void registerReaction(int rxnNumber, int type, int loc) { + m_index[rxnNumber] = std::pair(type, loc); + } + bool m_finalized; + }; +} + +#endif diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index c5fd96a2f..af8413426 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -191,78 +191,78 @@ namespace Cantera { } - /** - * Update properties that depend on concentrations. This method - * fills out the array of generalized concentrations by calling - * method getActivityConcentrations for each phase, which classes - * representing phases should overload to return the appropriate - * quantities. - */ - void InterfaceKinetics:: - _update_rates_C() { - int n; + /** + * Update properties that depend on concentrations. This method + * fills out the array of generalized concentrations by calling + * method getActivityConcentrations for each phase, which classes + * representing phases should overload to return the appropriate + * quantities. + */ + void InterfaceKinetics:: + _update_rates_C() { + int n; - int np = nPhases(); - for (n = 0; n < np; n++) { - /* - * We call the getActivityConcentrations function of each - * ThermoPhase class that makes up this kinetics object to - * obtain the generalized concentrations for species within that - * class. This is collected in the vector m_conc. m_start[] - * are integer indecises for that vector denoting the start of the - * species for each phase. - */ - thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]); - } - m_kdata->m_ROP_ok = false; + int np = nPhases(); + for (n = 0; n < np; n++) { + /* + * We call the getActivityConcentrations function of each + * ThermoPhase class that makes up this kinetics object to + * obtain the generalized concentrations for species within that + * class. This is collected in the vector m_conc. m_start[] + * are integer indecises for that vector denoting the start of the + * species for each phase. + */ + thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]); } + m_kdata->m_ROP_ok = false; + } - /** - * Update the equilibrium constants in molar units for all - * reversible reactions. Irreversible reactions have their - * equilibrium constant set to zero. - */ - void InterfaceKinetics::updateKc() { - int i, irxn; + /** + * Update the equilibrium constants in molar units for all + * reversible reactions. Irreversible reactions have their + * equilibrium constant set to zero. + */ + void InterfaceKinetics::updateKc() { + int i, irxn; - vector_fp& m_rkc = m_kdata->m_rkcn; - fill(m_rkc.begin(), m_rkc.end(), 0.0); + vector_fp& m_rkc = m_kdata->m_rkcn; + fill(m_rkc.begin(), m_rkc.end(), 0.0); - //static vector_fp mu(nTotalSpecies()); - if (m_nrev > 0) { + //static vector_fp mu(nTotalSpecies()); + if (m_nrev > 0) { - int n, nsp, k, ik=0; - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; - int np = nPhases(); - for (n = 0; n < np; n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); - nsp = thermo(n).nSpecies(); - for (k = 0; k < nsp; k++) { - m_mu0[ik] -= rt*thermo(n).logStandardConc(k); - m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); - ik++; - } - } + int n, nsp, k, ik = 0; + doublereal rt = GasConstant*thermo(0).temperature(); + doublereal rrt = 1.0 / rt; + int np = nPhases(); + for (n = 0; n < np; n++) { + thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); + nsp = thermo(n).nSpecies(); + for (k = 0; k < nsp; k++) { + m_mu0[ik] -= rt * thermo(n).logStandardConc(k); + m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); + ik++; + } + } - // compute Delta mu^0 for all reversible reactions - m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), - DATA_PTR(m_rkc)); + // compute Delta mu^0 for all reversible reactions + m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), + DATA_PTR(m_rkc)); - for (i = 0; i < m_nrev; i++) { - irxn = m_revindex[i]; - if (irxn < 0 || irxn >= nReactions()) { - throw CanteraError("InterfaceKinetics", - "illegal value: irxn = "+int2str(irxn)); - } - m_rkc[irxn] = exp(m_rkc[irxn]*rrt); - } - for (i = 0; i != m_nirrev; ++i) { - m_rkc[ m_irrev[i] ] = 0.0; - } - } + for (i = 0; i < m_nrev; i++) { + irxn = m_revindex[i]; + if (irxn < 0 || irxn >= nReactions()) { + throw CanteraError("InterfaceKinetics", + "illegal value: irxn = "+int2str(irxn)); + } + m_rkc[irxn] = exp(m_rkc[irxn]*rrt); + } + for (i = 0; i != m_nirrev; ++i) { + m_rkc[ m_irrev[i] ] = 0.0; + } } + } void InterfaceKinetics::checkPartialEquil() { diff --git a/Cantera/src/kinetics/KineticsFactory.cpp b/Cantera/src/kinetics/KineticsFactory.cpp index 770eef374..65d61d4db 100644 --- a/Cantera/src/kinetics/KineticsFactory.cpp +++ b/Cantera/src/kinetics/KineticsFactory.cpp @@ -22,6 +22,7 @@ #include "InterfaceKinetics.h" #include "EdgeKinetics.h" #include "importKinetics.h" +#include "AqueousKinetics.h" using namespace std; @@ -33,8 +34,8 @@ namespace Cantera { #endif static int ntypes = 5; - static string _types[] = {"none", "GasKinetics", "GRI30", "Interface", "Edge"}; - static int _itypes[] = {0, cGasKinetics, cGRI30, cInterfaceKinetics, cEdgeKinetics}; + static string _types[] = {"none", "GasKinetics", "GRI30", "Interface", "Edge", "AqueousKinetics"}; + static int _itypes[] = {0, cGasKinetics, cGRI30, cInterfaceKinetics, cEdgeKinetics, cAqueousKinetics}; /** * Return a new kinetics manager that implements a reaction @@ -105,6 +106,10 @@ namespace Cantera { case cEdgeKinetics: k = new EdgeKinetics; break; + + case cAqueousKinetics: + k = new AqueousKinetics; + break; default: throw UnknownKineticsModel("KineticsFactory::newKinetics", diff --git a/Cantera/src/kinetics/KineticsFactory.h b/Cantera/src/kinetics/KineticsFactory.h index 488b9255b..f4336c96d 100644 --- a/Cantera/src/kinetics/KineticsFactory.h +++ b/Cantera/src/kinetics/KineticsFactory.h @@ -31,7 +31,7 @@ namespace Cantera { CanteraError(proc, "Specified Kinetics model " + kineticsModel + " does not match any known type.") {} - virtual ~UnknownKineticsModel() {} + virtual ~UnknownKineticsModel() throw() {} }; diff --git a/Cantera/src/kinetics/Makefile.in b/Cantera/src/kinetics/Makefile.in index 843bbf93a..cfd625c43 100644 --- a/Cantera/src/kinetics/Makefile.in +++ b/Cantera/src/kinetics/Makefile.in @@ -36,13 +36,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) ifeq ($(do_kinetics),1) KINETICS_OBJ=importKinetics.o GRI_30_Kinetics.o KineticsFactory.o \ - GasKinetics.o \ + GasKinetics.o AqueousKinetics.o \ FalloffFactory.o ReactionStoichMgr.o Kinetics.o solveSP.o KINETICS_H = importKinetics.h GRI_30_Kinetics.h KineticsFactory.h \ Kinetics.h GasKinetics.h \ FalloffFactory.h ReactionStoichMgr.h reaction_defs.h \ FalloffMgr.h ThirdBodyMgr.h RateCoeffMgr.h ReactionData.h \ - RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h + RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h \ + AqueousKinetics.h KINETICS = $(KINETICS_OBJ) $(KINETICS_H) endif