From 8cbc08385454ddf714be7ce07e40b0c52e105c90 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 30 Mar 2012 23:48:20 +0000 Subject: [PATCH] Simplified AqueousKinetics by eliminating class AqueousKineticsData --- include/cantera/kinetics/AqueousKinetics.h | 63 +++------ src/kinetics/AqueousKinetics.cpp | 148 +++++++-------------- 2 files changed, 64 insertions(+), 147 deletions(-) diff --git a/include/cantera/kinetics/AqueousKinetics.h b/include/cantera/kinetics/AqueousKinetics.h index 104194872..f6aa11464 100644 --- a/include/cantera/kinetics/AqueousKinetics.h +++ b/include/cantera/kinetics/AqueousKinetics.h @@ -32,40 +32,7 @@ namespace Cantera // forward references - class ReactionData; -class AqueousKineticsData; -class Thermo; - -/** - * Holds mechanism-specific data. - */ -class AqueousKineticsData -{ -public: - AqueousKineticsData(); - - ~AqueousKineticsData(); - - AqueousKineticsData(const AqueousKineticsData& right); - - AqueousKineticsData& operator=(const AqueousKineticsData& right); - - doublereal m_logp_ref; - doublereal m_logc_ref; - vector_fp m_ropf; - vector_fp m_ropr; - vector_fp m_ropnet; - vector_fp m_rfn_low; - vector_fp m_rfn_high; - bool m_ROP_ok; - - doublereal m_temp; - vector_fp m_rfn; - - vector_fp m_rkcn; -}; - /** * Kinetics manager for elementary aqueous-phase chemistry. This @@ -141,7 +108,7 @@ public: */ virtual void getFwdRatesOfProgress(doublereal* fwdROP) { updateROP(); - std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP); + std::copy(m_ropf.begin(), m_ropf.end(), fwdROP); } /** @@ -152,7 +119,7 @@ public: */ virtual void getRevRatesOfProgress(doublereal* revROP) { updateROP(); - std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP); + std::copy(m_ropr.begin(), m_ropr.end(), revROP); } /** @@ -163,7 +130,7 @@ public: */ virtual void getNetRatesOfProgress(doublereal* netROP) { updateROP(); - std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); + std::copy(m_ropnet.begin(), m_ropnet.end(), netROP); } @@ -253,8 +220,7 @@ public: //#ifdef HWMECH //get_wdot(&m_kdata->m_ropnet[0], net); //#else - m_rxnstoich->getNetProductionRates(m_kk, - &m_kdata->m_ropnet[0], net); + m_rxnstoich->getNetProductionRates(m_kk, &m_ropnet[0], net); //#endif } @@ -267,8 +233,7 @@ public: */ virtual void getCreationRates(doublereal* cdot) { updateROP(); - m_rxnstoich->getCreationRates(m_kk, &m_kdata->m_ropf[0], - &m_kdata->m_ropr[0], cdot); + m_rxnstoich->getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot); } /** @@ -280,9 +245,7 @@ public: */ virtual void getDestructionRates(doublereal* ddot) { updateROP(); - m_rxnstoich->getDestructionRates(m_kk, &m_kdata->m_ropf[0], - &m_kdata->m_ropr[0], ddot); - + m_rxnstoich->getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot); } //@} @@ -406,11 +369,21 @@ protected: std::vector m_rxneqn; - AqueousKineticsData* m_kdata; - vector_fp m_conc; vector_fp m_grt; + //! @name Aqueous kinetics data + //!@{ + vector_fp m_ropf; + vector_fp m_ropr; + vector_fp m_ropnet; + bool m_ROP_ok; + + doublereal m_temp; + vector_fp m_rfn; + + vector_fp m_rkcn; + //!@} private: diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index 32bca1a1c..4c31a71fd 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -21,49 +21,6 @@ using namespace std; namespace Cantera { -//==================================================================================================================== -AqueousKineticsData::AqueousKineticsData() : - m_logp_ref(0.0), - m_logc_ref(0.0), - m_ROP_ok(false), - m_temp(0.0) -{ -} -//==================================================================================================================== -AqueousKineticsData::~AqueousKineticsData() -{ -} -//==================================================================================================================== -AqueousKineticsData::AqueousKineticsData(const AqueousKineticsData& right) : - m_logp_ref(0.0), - m_logc_ref(0.0), - m_ROP_ok(false), - m_temp(0.0) -{ - *this=right; -} -//==================================================================================================================== - -AqueousKineticsData& AqueousKineticsData::operator=(const AqueousKineticsData& right) -{ - if (this != &right) { - m_logp_ref = right.m_logp_ref; - m_logc_ref = right.m_logc_ref; - m_ropf = right.m_ropf; - m_ropr = right.m_ropr; - m_ropnet = right.m_ropnet; - m_rfn_low = right.m_rfn_low; - m_rfn_high = right.m_rfn_high; - m_ROP_ok = right.m_ROP_ok; - m_temp = right.m_temp; - m_rfn = right.m_rfn; - m_rkcn = right.m_rkcn; - } - return *this; -} -//==================================================================================================================== - - //==================================================================================================================== /** * Construct an empty reaction mechanism. @@ -73,13 +30,13 @@ AqueousKinetics::AqueousKinetics(thermo_t* thermo) : m_nfall(0), m_nirrev(0), m_nrev(0), + m_ROP_ok(false), + m_temp(0.0), m_finalized(false) { if (thermo != 0) { addPhase(*thermo); } - m_kdata = new AqueousKineticsData; - m_kdata->m_temp = 0.0; m_rxnstoich = new ReactionStoichMgr; } //==================================================================================================================== @@ -88,6 +45,8 @@ AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) : m_nfall(0), m_nirrev(0), m_nrev(0), + m_ROP_ok(false), + m_temp(0.0), m_finalized(false) { *this = right; @@ -95,7 +54,6 @@ AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) : //==================================================================================================================== AqueousKinetics::~AqueousKinetics() { - delete m_kdata; delete m_rxnstoich; } //==================================================================================================================== @@ -126,7 +84,13 @@ AqueousKinetics& AqueousKinetics::operator=(const AqueousKinetics& right) m_revindex = right.m_revindex; m_rxneqn = right.m_rxneqn; - *m_kdata = *(right.m_kdata); + m_ropf = right.m_ropf; + m_ropr = right.m_ropr; + m_ropnet = right.m_ropnet; + m_ROP_ok = right.m_ROP_ok; + m_temp = right.m_temp; + m_rfn = right.m_rfn; + m_rkcn = right.m_rkcn; m_conc = right.m_conc; m_grt = right.m_grt; @@ -160,15 +124,12 @@ 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_rates.update(T, logT, &m_rfn[0]); - - m_kdata->m_temp = T; + m_temp = T; updateKc(); - m_kdata->m_ROP_ok = false; - + m_ROP_ok = false; }; @@ -181,7 +142,7 @@ _update_rates_C() { thermo().getActivityConcentrations(&m_conc[0]); - m_kdata->m_ROP_ok = false; + m_ROP_ok = false; } /** @@ -189,28 +150,27 @@ _update_rates_C() */ void AqueousKinetics::updateKc() { - vector_fp& m_rkc = m_kdata->m_rkcn; - doublereal rt = GasConstant* m_kdata->m_temp; + doublereal rt = GasConstant * m_temp; thermo().getStandardChemPotentials(&m_grt[0]); - fill(m_rkc.begin(), m_rkc.end(), 0.0); + fill(m_rkcn.begin(), m_rkcn.end(), 0.0); for (size_t k = 0; k < thermo().nSpecies(); 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]); + m_rxnstoich->getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]); //doublereal logStandConc = m_kdata->m_logStandConc; doublereal rrt = 1.0/(GasConstant * thermo().temperature()); for (size_t i = 0; i < m_nrev; i++) { size_t irxn = m_revindex[i]; - m_rkc[irxn] = exp(m_rkc[irxn]*rrt); + m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt); } for (size_t i = 0; i != m_nirrev; ++i) { - m_rkc[ m_irrev[i] ] = 0.0; + m_rkcn[ m_irrev[i] ] = 0.0; } } @@ -221,27 +181,26 @@ void AqueousKinetics::updateKc() void AqueousKinetics::getEquilibriumConstants(doublereal* kc) { _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; + fill(m_rkcn.begin(), m_rkcn.end(), 0.0); + doublereal rt = GasConstant * m_temp; for (size_t k = 0; k < thermo().nSpecies(); 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]); + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]); doublereal rrt = 1.0/(GasConstant * thermo().temperature()); for (size_t i = 0; i < m_ii; i++) { - kc[i] = exp(-rkc[i]*rrt); + kc[i] = exp(-m_rkcn[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; + m_temp = 0.0; } /** @@ -411,49 +370,41 @@ void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS) void AqueousKinetics::updateROP() { - _update_rates_T(); _update_rates_C(); - if (m_kdata->m_ROP_ok) { + if (m_ROP_ok) { return; } - const vector_fp& rf = m_kdata->m_rfn; - const vector_fp& m_rkc = m_kdata->m_rkcn; - vector_fp& ropf = m_kdata->m_ropf; - vector_fp& ropr = m_kdata->m_ropr; - vector_fp& ropnet = m_kdata->m_ropnet; - // copy rate coefficients into ropf - copy(rf.begin(), rf.end(), ropf.begin()); - + copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); // multiply by perturbation factor - multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); + multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin()); // copy the forward rates to the reverse rates - copy(ropf.begin(), ropf.end(), ropr.begin()); + 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(ropr.begin(), ropr.end(), m_rkc.begin()); + multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); // multiply ropf by concentration products - m_rxnstoich->multiplyReactants(&m_conc[0], &ropf[0]); + m_rxnstoich->multiplyReactants(&m_conc[0], &m_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_rxnstoich->multiplyRevProducts(&m_conc[0], &m_ropr[0]); //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); for (size_t j = 0; j != m_ii; ++j) { - ropnet[j] = ropf[j] - ropr[j]; + m_ropnet[j] = m_ropf[j] - m_ropr[j]; } - m_kdata->m_ROP_ok = true; + m_ROP_ok = true; } /** @@ -471,17 +422,13 @@ getFwdRateConstants(doublereal* kfwd) _update_rates_C(); // copy rate coefficients into ropf - const vector_fp& rf = m_kdata->m_rfn; - vector_fp& ropf = m_kdata->m_ropf; - copy(rf.begin(), rf.end(), ropf.begin()); - - + copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); // multiply by perturbation factor - multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); + multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin()); for (size_t i = 0; i < m_ii; i++) { - kfwd[i] = ropf[i]; + kfwd[i] = m_ropf[i]; } } @@ -507,18 +454,16 @@ getRevRateConstants(doublereal* krev, bool doIrreversible) getFwdRateConstants(krev); if (doIrreversible) { - doublereal* tmpKc = &m_kdata->m_ropnet[0]; - getEquilibriumConstants(tmpKc); + getEquilibriumConstants(&m_ropnet[0]); for (size_t i = 0; i < m_ii; i++) { - krev[i] /= tmpKc[i]; + krev[i] /= m_ropnet[i]; } } else { /* - * m_rkc[] is zero for irreversibly reactions + * m_rkcn[] is zero for irreversible reactions */ - const vector_fp& m_rkc = m_kdata->m_rkcn; for (size_t i = 0; i < m_ii; i++) { - krev[i] *= m_rkc[i]; + krev[i] *= m_rkcn[i]; } } } @@ -548,7 +493,7 @@ void AqueousKinetics::addElementaryReaction(ReactionData& r) iloc = m_rates.install(reactionNumber(), r); // add constant term to rate coeff value vector - m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); + m_rfn.push_back(r.rateCoeffParameters[0]); // forward rxn order equals number of reactants m_fwdOrder.push_back(r.reactants.size()); @@ -561,9 +506,9 @@ void AqueousKinetics::addElementaryReaction(ReactionData& r) 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); + m_ropf.push_back(0.0); // extend by one for new rxn + m_ropr.push_back(0.0); + m_ropnet.push_back(0.0); size_t n, ns, m; doublereal nsFlt; doublereal reactantGlobalOrder = 0.0; @@ -610,7 +555,7 @@ void AqueousKinetics::installReagents(const ReactionData& r) } m_products.push_back(pk); - m_kdata->m_rkcn.push_back(0.0); + m_rkcn.push_back(0.0); m_rxnstoich->add(reactionNumber(), r); @@ -645,7 +590,6 @@ void AqueousKinetics::init() 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()