From 7d901ed74f1255a33f26c0c985b7fc0b5f25e71c Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 25 Jun 2010 16:04:29 +0000 Subject: [PATCH] Added in duplication routines. Added in a shallow pointer copy facility. --- Cantera/src/kinetics/GasKinetics.cpp | 1175 +++++++++++--------- Cantera/src/kinetics/GasKinetics.h | 748 +++++++------ Cantera/src/kinetics/InterfaceKinetics.cpp | 10 +- Cantera/src/kinetics/InterfaceKinetics.h | 14 +- Cantera/src/kinetics/Kinetics.cpp | 72 +- Cantera/src/kinetics/Kinetics.h | 60 +- 6 files changed, 1134 insertions(+), 945 deletions(-) diff --git a/Cantera/src/kinetics/GasKinetics.cpp b/Cantera/src/kinetics/GasKinetics.cpp index 4633a08ee..9aeff996e 100644 --- a/Cantera/src/kinetics/GasKinetics.cpp +++ b/Cantera/src/kinetics/GasKinetics.cpp @@ -29,636 +29,711 @@ using namespace std; namespace Cantera { + //==================================================================================================================== + /* + * Construct an empty reaction mechanism. + */ + GasKinetics:: + GasKinetics(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 GasKineticsData; + m_kdata->m_temp = 0.0; + m_rxnstoich = new ReactionStoichMgr; + } - /** - * Construct an empty reaction mechanism. - */ - GasKinetics:: - GasKinetics(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 GasKineticsData; - m_kdata->m_temp = 0.0; - m_rxnstoich = new ReactionStoichMgr; - } + //==================================================================================================================== + GasKinetics::GasKinetics(const GasKinetics &right) : + Kinetics(), + m_kk(0), + m_nfall(0), + m_nirrev(0), + m_nrev(0), + m_finalized(false) + { + *this = right; + } + //==================================================================================================================== + GasKinetics::~GasKinetics() + { + delete m_kdata; + delete m_rxnstoich; + } + //==================================================================================================================== + GasKinetics& GasKinetics::operator=(const GasKinetics &right) + { + if (this == &right) return *this; - GasKinetics:: - ~GasKinetics() {delete m_kdata; delete m_rxnstoich;} + Kinetics::operator=(right); + + m_kk = right.m_kk; + m_nfall = right.m_nfall; + m_fallindx = right.m_fallindx; + m_falloff_low_rates = right.m_falloff_low_rates; + m_falloff_high_rates = right.m_falloff_high_rates; + m_rates = right.m_rates; + m_index = right.m_index; + m_falloffn = right.m_falloffn; + m_3b_concm = right.m_3b_concm; + m_falloff_concm = right.m_falloff_concm; + m_irrev = right.m_irrev; + m_rxnstoich = right.m_rxnstoich; + m_fwdOrder = right.m_fwdOrder; + m_nirrev = right.m_nirrev; + m_nrev = right.m_nrev; + m_rgroups = right.m_rgroups; + m_pgroups = right.m_pgroups; + m_rxntype = right.m_rxntype; + m_rrxn = right.m_rrxn; + m_prxn = right.m_prxn; + m_dn = right.m_dn; + m_revindex = right.m_revindex; + m_rxneqn = right.m_rxneqn; + m_kdata = right.m_kdata; + m_conc = right.m_conc; + m_grt = right.m_grt; + m_finalized = right.m_finalized; - /** - * Update temperature-dependent portions of reaction rates and - * falloff functions. - */ - void GasKinetics:: - update_T() {} + throw CanteraError("GasKinetics::operator=()", + "Unfinished implementation"); - void GasKinetics:: - update_C() {} + return *this; + } + //==================================================================================================================== + // Duplication routine for objects which inherit from Kinetics + /* + * This virtual routine can be used to duplicate %Kinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object + */ + Kinetics *GasKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const { + GasKinetics* gK = new GasKinetics(*this); + gK->assignShallowPointers(tpVector); + return dynamic_cast(gK); + } + //==================================================================================================================== + /** + * Update temperature-dependent portions of reaction rates and + * falloff functions. + */ + void GasKinetics::update_T() + { + } + //==================================================================================================================== + void GasKinetics:: + update_C() {} + //==================================================================================================================== + void GasKinetics:: + _update_rates_T() { + doublereal T = thermo().temperature(); + m_kdata->m_logStandConc = log(thermo().standardConcentration()); + //if (fabs(T - m_kdata->m_temp) > 0.0) { + doublereal logT = log(T); + m_rates.update(T, logT, &m_kdata->m_rfn[0]); + m_falloff_low_rates.update(T, logT, &m_kdata->m_rfn_low[0]); + m_falloff_high_rates.update(T, logT, &m_kdata->m_rfn_high[0]); + m_falloffn.updateTemp(T, &m_kdata->falloff_work[0]); + m_kdata->m_temp = T; + updateKc(); + m_kdata->m_ROP_ok = false; + //} + }; - void GasKinetics:: - _update_rates_T() { - doublereal T = thermo().temperature(); - m_kdata->m_logStandConc = log(thermo().standardConcentration()); - //if (fabs(T - m_kdata->m_temp) > 0.0) { - doublereal logT = log(T); - m_rates.update(T, logT, &m_kdata->m_rfn[0]); - m_falloff_low_rates.update(T, logT, &m_kdata->m_rfn_low[0]); - m_falloff_high_rates.update(T, logT, &m_kdata->m_rfn_high[0]); - m_falloffn.updateTemp(T, &m_kdata->falloff_work[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 GasKinetics:: - _update_rates_C() { - thermo().getActivityConcentrations(&m_conc[0]); - doublereal ctot = thermo().molarDensity(); - m_3b_concm.update(m_conc, ctot, &m_kdata->concm_3b_values[0]); - m_falloff_concm.update(m_conc, ctot, - &m_kdata->concm_falloff_values[0]); - m_kdata->m_ROP_ok = false; - } - - /** - * Update the equilibrium constants in molar units. - */ - void GasKinetics::updateKc() { - int i, irxn; - vector_fp& m_rkc = m_kdata->m_rkcn; + //==================================================================================================================== + /** + * Update properties that depend on concentrations. Currently only + * the enhanced collision partner concentrations are updated here. + */ + void GasKinetics:: + _update_rates_C() { + thermo().getActivityConcentrations(&m_conc[0]); + doublereal ctot = thermo().molarDensity(); + m_3b_concm.update(m_conc, ctot, &m_kdata->concm_3b_values[0]); + m_falloff_concm.update(m_conc, ctot, + &m_kdata->concm_falloff_values[0]); + m_kdata->m_ROP_ok = false; + } + //==================================================================================================================== + /** + * Update the equilibrium constants in molar units. + */ + void GasKinetics::updateKc() { + int i, irxn; + vector_fp& m_rkc = m_kdata->m_rkcn; - thermo().getStandardChemPotentials(&m_grt[0]); - fill(m_rkc.begin(), m_rkc.end(), 0.0); + thermo().getStandardChemPotentials(&m_grt[0]); + fill(m_rkc.begin(), m_rkc.end(), 0.0); - // compute Delta G^0 for all reversible reactions - m_rxnstoich->getRevReactionDelta(m_ii, &m_grt[0], &m_rkc[0]); + // 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 - m_dn[irxn]*logStandConc); - } - - for(i = 0; i != m_nirrev; ++i) { - m_rkc[ m_irrev[i] ] = 0.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 - m_dn[irxn]*logStandConc); } - /** - * Get the equilibrium constants of all reactions, whether - * reversible or not. - */ - void GasKinetics::getEquilibriumConstants(doublereal* kc) { - int i; - _update_rates_T(); - vector_fp& rkc = m_kdata->m_rkcn; - //thermo().getGibbs_RT(m_grt.begin()); - thermo().getStandardChemPotentials(&m_grt[0]); - fill(rkc.begin(), rkc.end(), 0.0); + 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 GasKinetics::getEquilibriumConstants(doublereal* kc) { + int i; + _update_rates_T(); + vector_fp& rkc = m_kdata->m_rkcn; + //thermo().getGibbs_RT(m_grt.begin()); + thermo().getStandardChemPotentials(&m_grt[0]); + fill(rkc.begin(), rkc.end(), 0.0); - // compute Delta G^0 for all reactions - m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], &rkc[0]); + // compute Delta G^0 for all reactions + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], &rkc[0]); - doublereal logStandConc = m_kdata->m_logStandConc; - doublereal rrt = 1.0/(GasConstant * thermo().temperature()); - for (i = 0; i < m_ii; i++) { - kc[i] = exp(-rkc[i]*rrt + m_dn[i]*logStandConc); - } - - // 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; + doublereal logStandConc = m_kdata->m_logStandConc; + doublereal rrt = 1.0/(GasConstant * thermo().temperature()); + for (i = 0; i < m_ii; i++) { + kc[i] = exp(-rkc[i]*rrt + m_dn[i]*logStandConc); } - /** - * - * 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 + // 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 GasKinetics::getDeltaGibbs(doublereal* deltaG) { + /* + * Get the chemical potentials of the species in the + * ideal gas solution. */ - void GasKinetics::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 + thermo().getChemPotentials(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. */ - void GasKinetics::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 + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaG); + } + //==================================================================================================================== + /** + * + * getDeltaEnthalpy(): + * + * Return the vector of values for the reactions change in + * enthalpy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaEnthalpy(doublereal* deltaH) { + /* + * Get the partial molar enthalpy of all species in the + * ideal gas. */ - void GasKinetics::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 + thermo().getPartialMolarEnthalpies(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. */ - void GasKinetics::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 + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaH); + } + //==================================================================================================================== + /* + * + * getDeltaEntropy(): + * + * Return the vector of values for the reactions change in + * entropy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + void GasKinetics::getDeltaEntropy( doublereal* deltaS) { + /* + * Get the partial molar entropy of all species in the + * solid solution. */ - void GasKinetics::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 + thermo().getPartialMolarEntropies(&m_grt[0]); + /* + * Use the stoichiometric manager to find deltaS for each + * reaction. */ - void GasKinetics::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); + m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaS); + } + //==================================================================================================================== + /** + * + * getDeltaSSGibbs(): + * + * Return the vector of values for the reaction + * standard state gibbs free energy change. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaSSGibbs(doublereal* deltaG) { + /* + * 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 GasKinetics::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 GasKinetics::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 GasKinetics::processFalloffReactions() { + + int i; + const vector_fp& fc = m_kdata->concm_falloff_values; + const array_fp& m_rf_low = m_kdata->m_rfn_low; + const array_fp& m_rf_high = m_kdata->m_rfn_high; + + // use m_ropr for temporary storage of reduced pressure + array_fp& pr = m_kdata->m_ropr; + + array_fp& ropf = m_kdata->m_ropf; + + for (i = 0; i < m_nfall; i++) { + pr[i] = fc[i] * m_rf_low[i] / m_rf_high[i]; } - void GasKinetics::processFalloffReactions() { - - int i; - const vector_fp& fc = m_kdata->concm_falloff_values; - const array_fp& m_rf_low = m_kdata->m_rfn_low; - const array_fp& m_rf_high = m_kdata->m_rfn_high; - - // use m_ropr for temporary storage of reduced pressure - array_fp& pr = m_kdata->m_ropr; - - array_fp& ropf = m_kdata->m_ropf; - - for (i = 0; i < m_nfall; i++) { - pr[i] = fc[i] * m_rf_low[i] / m_rf_high[i]; - } - - m_falloffn.pr_to_falloff( &pr[0], &m_kdata->falloff_work[0] ); + m_falloffn.pr_to_falloff( &pr[0], &m_kdata->falloff_work[0] ); - for (i = 0; i < m_nfall; i++) { - pr[i] *= m_rf_high[i]; - } - - scatter_copy(pr.begin(), pr.begin() + m_nfall, - ropf.begin(), m_fallindx.begin()); + for (i = 0; i < m_nfall; i++) { + pr[i] *= m_rf_high[i]; } + scatter_copy(pr.begin(), pr.begin() + m_nfall, + ropf.begin(), m_fallindx.begin()); + } - void GasKinetics::updateROP() { + //==================================================================================================================== + void GasKinetics::updateROP() { - _update_rates_T(); - _update_rates_C(); + _update_rates_T(); + _update_rates_C(); - if (m_kdata->m_ROP_ok) return; + 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; + 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()); + // copy rate coefficients into ropf + copy(rf.begin(), rf.end(), ropf.begin()); - // multiply ropf by enhanced 3b conc for all 3b rxns - m_3b_concm.multiply( &ropf[0], &m_kdata->concm_3b_values[0] ); + // multiply ropf by enhanced 3b conc for all 3b rxns + m_3b_concm.multiply( &ropf[0], &m_kdata->concm_3b_values[0] ); - processFalloffReactions(); + processFalloffReactions(); - // multiply by perturbation factor - multiply_each(ropf.begin(), ropf.end(), m_perturb.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()); + // 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()); + // 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()); + // 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 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; + for (int j = 0; j != m_ii; ++j) { + ropnet[j] = ropf[j] - ropr[j]; } - /** - * - * 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. + 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 GasKinetics:: + 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 ropf by enhanced 3b conc for all 3b rxns + m_3b_concm.multiply(&ropf[0], &m_kdata->concm_3b_values[0] ); + + /* + * This routine is hardcoded to replace some of the values + * of the ropf vector. */ - void GasKinetics:: - getFwdRateConstants(doublereal *kfwd) { - _update_rates_T(); - _update_rates_C(); + processFalloffReactions(); - // 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 ropf by enhanced 3b conc for all 3b rxns - m_3b_concm.multiply(&ropf[0], &m_kdata->concm_3b_values[0] ); - - /* - * This routine is hardcoded to replace some of the values - * of the ropf vector. - */ - processFalloffReactions(); - - // multiply by perturbation factor - multiply_each(ropf.begin(), ropf.end(), m_perturb.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]; - } + 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. + } + //==================================================================================================================== + /** + * + * getRevRateConstants(): + * + * Return a vector of the reverse reaction rate constants + * + * Length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. + */ + void GasKinetics:: + getRevRateConstants(doublereal *krev, bool doIrreversible) { + /* + * go get the forward rate constants. -> note, we don't + * really care about speed or redundancy in these + * informational routines. */ - void GasKinetics:: - 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); + 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]; - } - } + 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 GasKinetics:: + addReaction(const ReactionData& r) { - void GasKinetics:: - addReaction(const ReactionData& r) { + if (r.reactionType == ELEMENTARY_RXN) addElementaryReaction(r); + else if (r.reactionType == THREE_BODY_RXN) addThreeBodyReaction(r); + else if (r.reactionType == FALLOFF_RXN) addFalloffReaction(r); - if (r.reactionType == ELEMENTARY_RXN) addElementaryReaction(r); - else if (r.reactionType == THREE_BODY_RXN) addThreeBodyReaction(r); - else if (r.reactionType == FALLOFF_RXN) addFalloffReaction(r); + // operations common to all reaction types + installReagents( r ); + installGroups(reactionNumber(), r.rgroups, r.pgroups); + incrementRxnCount(); + m_rxneqn.push_back(r.equation); + } - // operations common to all reaction types - installReagents( r ); - installGroups(reactionNumber(), r.rgroups, r.pgroups); - incrementRxnCount(); - m_rxneqn.push_back(r.equation); - } + //==================================================================================================================== + void GasKinetics:: + addFalloffReaction(const ReactionData& r) { + // install high and low rate coeff calculators - void GasKinetics:: - addFalloffReaction(const ReactionData& r) { - - // install high and low rate coeff calculators - - int iloc = m_falloff_high_rates.install(m_nfall, - r.rateCoeffType, - r.rateCoeffParameters.size(), - &r.rateCoeffParameters[0] ); + int iloc = m_falloff_high_rates.install(m_nfall, + r.rateCoeffType, + r.rateCoeffParameters.size(), + &r.rateCoeffParameters[0] ); - m_falloff_low_rates.install( m_nfall, - r.rateCoeffType, r.auxRateCoeffParameters.size(), - DATA_PTR(r.auxRateCoeffParameters) ); + m_falloff_low_rates.install( m_nfall, + r.rateCoeffType, r.auxRateCoeffParameters.size(), + DATA_PTR(r.auxRateCoeffParameters) ); - // add constant terms to high and low rate - // coeff value vectors - m_kdata->m_rfn_high.push_back(r.rateCoeffParameters[0]); - m_kdata->m_rfn_low.push_back(r.auxRateCoeffParameters[0]); + // add constant terms to high and low rate + // coeff value vectors + m_kdata->m_rfn_high.push_back(r.rateCoeffParameters[0]); + m_kdata->m_rfn_low.push_back(r.auxRateCoeffParameters[0]); - // add a dummy entry in m_rf, where computed falloff - // rate coeff will be put - m_kdata->m_rfn.push_back(0.0); + // add a dummy entry in m_rf, where computed falloff + // rate coeff will be put + m_kdata->m_rfn.push_back(0.0); - // add this reaction number to the list of - // falloff reactions - m_fallindx.push_back( reactionNumber() ); + // add this reaction number to the list of + // falloff reactions + m_fallindx.push_back( reactionNumber() ); - // install the enhanced third-body concentration - // calculator for this reaction - m_falloff_concm.install( m_nfall, r.thirdBodyEfficiencies, - r.default_3b_eff); + // install the enhanced third-body concentration + // calculator for this reaction + m_falloff_concm.install( m_nfall, r.thirdBodyEfficiencies, + r.default_3b_eff); - // install the falloff function calculator for - // this reaction - m_falloffn.install( m_nfall, r.falloffType, r.falloffParameters ); + // install the falloff function calculator for + // this reaction + m_falloffn.install( m_nfall, r.falloffType, r.falloffParameters ); - // forward rxn order equals number of reactants, since rate - // coeff is defined in terms of the high-pressure limit - m_fwdOrder.push_back(r.reactants.size()); + // forward rxn order equals number of reactants, since rate + // coeff is defined in terms of the high-pressure limit + m_fwdOrder.push_back(r.reactants.size()); - // increment the falloff reaction counter - ++m_nfall; - registerReaction( reactionNumber(), FALLOFF_RXN, iloc); - } + // increment the falloff reaction counter + ++m_nfall; + registerReaction( reactionNumber(), FALLOFF_RXN, iloc); + } + //==================================================================================================================== + void GasKinetics:: + addElementaryReaction(const ReactionData& r) { + int iloc; - void GasKinetics:: - addElementaryReaction(const ReactionData& r) { - int iloc; + // install rate coeff calculator + iloc = m_rates.install( reactionNumber(), + r.rateCoeffType, r.rateCoeffParameters.size(), + DATA_PTR(r.rateCoeffParameters) ); - // 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]); - // 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); + } - // forward rxn order equals number of reactants - m_fwdOrder.push_back(r.reactants.size()); - registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc); - } - - - void GasKinetics:: - addThreeBodyReaction(const ReactionData& r) { + //==================================================================================================================== + void GasKinetics:: + addThreeBodyReaction(const ReactionData& r) { - int iloc; - // install rate coeff calculator - iloc = m_rates.install( reactionNumber(), - r.rateCoeffType, r.rateCoeffParameters.size(), - DATA_PTR(r.rateCoeffParameters) ); + 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]); + // add constant term to rate coeff value vector + m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); - // forward rxn order equals number of reactants + 1 - m_fwdOrder.push_back(r.reactants.size() + 1); - - m_3b_concm.install( reactionNumber(), r.thirdBodyEfficiencies, - r.default_3b_eff ); - registerReaction( reactionNumber(), THREE_BODY_RXN, iloc); - } + // forward rxn order equals number of reactants + 1 + m_fwdOrder.push_back(r.reactants.size() + 1); + m_3b_concm.install( reactionNumber(), r.thirdBodyEfficiencies, + r.default_3b_eff ); + registerReaction( reactionNumber(), THREE_BODY_RXN, iloc); + } + //==================================================================================================================== - void GasKinetics::installReagents(const ReactionData& r) { + void GasKinetics::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(); + 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++; - } + 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); - - void GasKinetics::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; - } + 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); - void GasKinetics::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); + 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 GasKinetics::finalize() { - if (!m_finalized) { - // int i, j, nr, np; - m_kdata->falloff_work.resize( - static_cast(m_falloffn.workSize())); - m_kdata->concm_3b_values.resize( - static_cast(m_3b_concm.workSize())); - m_kdata->concm_falloff_values.resize( - static_cast(m_falloff_concm.workSize())); - -// for (i = 0; i < m_ii; i++) { -// nr = m_reactants[i].size(); -// for (j = 0; j < nr; j++) { -// m_rstoich[i][m_reactants[i][j]]++; -// } -// np = m_products[i].size(); -// for (j = 0; j < np; j++) { -// m_pstoich[i][m_products[i][j]]++; -// } -// } - //m_rxnstoich->write("c.cpp"); - m_finalized = true; - } + void GasKinetics::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; } + } - bool GasKinetics::ready() const { - return (m_finalized); + //==================================================================================================================== + void GasKinetics::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 GasKinetics::finalize() { + if (!m_finalized) { + // int i, j, nr, np; + m_kdata->falloff_work.resize( + static_cast(m_falloffn.workSize())); + m_kdata->concm_3b_values.resize( + static_cast(m_3b_concm.workSize())); + m_kdata->concm_falloff_values.resize( + static_cast(m_falloff_concm.workSize())); + + // for (i = 0; i < m_ii; i++) { + // nr = m_reactants[i].size(); + // for (j = 0; j < nr; j++) { + // m_rstoich[i][m_reactants[i][j]]++; + // } + // np = m_products[i].size(); + // for (j = 0; j < np; j++) { + // m_pstoich[i][m_products[i][j]]++; + // } + // } + //m_rxnstoich->write("c.cpp"); + m_finalized = true; } - + } + //==================================================================================================================== + bool GasKinetics::ready() const { + return (m_finalized); + } + //==================================================================================================================== } +//====================================================================================================================== diff --git a/Cantera/src/kinetics/GasKinetics.h b/Cantera/src/kinetics/GasKinetics.h index 7dee80dcd..f221a76ab 100644 --- a/Cantera/src/kinetics/GasKinetics.h +++ b/Cantera/src/kinetics/GasKinetics.h @@ -34,385 +34,435 @@ void get_wdot(const doublereal* rop, doublereal* wdot); namespace Cantera { - // forward references + // forward references - class Enhanced3BConc; - class ReactionData; - class GasKineticsData; - class Thermo; + class Enhanced3BConc; + class ReactionData; + class GasKineticsData; + class Thermo; + + /** + * Holds mechanism-specific data. + */ + class GasKineticsData { + public: + GasKineticsData() : + m_logp_ref(0.0), + m_logc_ref(0.0), + m_logStandConc(0.0), + m_ROP_ok(false), + m_temp(0.0) + {} + virtual ~GasKineticsData(){} + + doublereal m_logp_ref, m_logc_ref, m_logStandConc; + array_fp m_ropf, 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 falloff_work; + array_fp concm_3b_values; + array_fp concm_falloff_values; + array_fp m_rkcn; + }; + + + /** + * Kinetics manager for elementary gas-phase chemistry. This + * kinetics manager implements standard mass-action reaction rate + * expressions for low-density gases. + * @ingroup kinetics + */ + class GasKinetics : public Kinetics { + + public: /** - * Holds mechanism-specific data. + * @name Constructors and General Information */ - class GasKineticsData { - public: - GasKineticsData() : - m_logp_ref(0.0), - m_logc_ref(0.0), - m_logStandConc(0.0), - m_ROP_ok(false), - m_temp(0.0) - {} - virtual ~GasKineticsData(){} + //@{ - doublereal m_logp_ref, m_logc_ref, m_logStandConc; - array_fp m_ropf, m_ropr, m_ropnet; - array_fp m_rfn_low, m_rfn_high; - bool m_ROP_ok; + //! Constructor. + /*! + * @param thermo Pointer to the gas ThermoPhase (optional) + */ + GasKinetics(thermo_t* thermo = 0); - doublereal m_temp; - array_fp m_rfn; - array_fp falloff_work; - array_fp concm_3b_values; - array_fp concm_falloff_values; - array_fp m_rkcn; - }; + + //!Copy Constructor for the %GasKinetics object. + /*! + * Currently, this is not fully implemented. If called it will + * throw an exception. + * + * @param right object to be copied + */ + GasKinetics(const GasKinetics &right); + + //! Destructor. + virtual ~GasKinetics(); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %GasKinetics object to be copied into the + * current one. + */ + GasKinetics& operator=(const GasKinetics &right); + + //! Duplication routine for objects which inherit from Kinetics + /*! + * This virtual routine can be used to duplicate %Kinetics objects + * inherited from %Kinetics even if the application only has + * a pointer to %Kinetics to work with. + * + * These routines are basically wrappers around the derived copy constructor. + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object + */ + virtual Kinetics *duplMyselfAsKinetics(const std::vector & tpVector) const; + + + //! Identifies the subclass of the Kinetics manager type. + /*! + * These are listed in mix_defs.h. + */ + virtual int ID() const { return cGasKinetics; } + + //! Identifies the kinetics manager type. + /*! + * Each class derived from Kinetics should overload this method to + * return a unique integer. Standard values are defined in file + * mix_defs.h. + */ + virtual int type() const { return cGasKinetics; } + + 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); + } /** - * Kinetics manager for elementary gas-phase chemistry. This - * kinetics manager implements standard mass-action reaction rate - * expressions for low-density gases. - * @ingroup kinetics + * 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. */ - class GasKinetics : public Kinetics { + virtual void getEquilibriumConstants(doublereal* kc); - public: + /** + * 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); - /** - * @name Constructors and General Information - */ - //@{ - /// Constructor. - GasKinetics(thermo_t* thermo = 0); + /** + * 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); - /// Destructor. - virtual ~GasKinetics(); + /** + * 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); - virtual int ID() const { return cGasKinetics; } - virtual int type() const { return cGasKinetics; } + /** + * 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); - virtual doublereal reactantStoichCoeff(int k, int i) const { - return m_rrxn[k][i]; - } + /** + * 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); - virtual doublereal productStoichCoeff(int k, int i) const { - return m_prxn[k][i]; - } + /** + * 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 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); - } + //@} + /** + * @name Species Production Rates + */ + //@{ - /** - * 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); - } + //! 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 + } - /** - * 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); - } + /** + * 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); + // fill(ddot, ddot + m_kk, 0.0); + //m_revProductStoich.incrementSpecies( + // m_kdata->m_ropr.begin(), ddot); + //m_reactantStoich.incrementSpecies( + // m_kdata->m_ropf.begin(), 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(); - /** - * 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); - // fill(ddot, ddot + m_kk, 0.0); - //m_revProductStoich.incrementSpecies( - // m_kdata->m_ropr.begin(), ddot); - //m_reactantStoich.incrementSpecies( - // m_kdata->m_ropf.begin(), 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]; } - 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; + + int m_nfall; + + array_int m_fallindx; + + Rate1 m_falloff_low_rates; + Rate1 m_falloff_high_rates; + Rate1 m_rates; + + mutable std::map > m_index; + + FalloffMgr m_falloffn; + + ThirdBodyMgr m_3b_concm; + ThirdBodyMgr m_falloff_concm; + + 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; + + GasKineticsData* m_kdata; + + array_fp m_conc; + void processFalloffReactions(); + array_fp m_grt; - void _update_rates_T(); - void _update_rates_C(); + private: - //@} + int reactionNumber(){ return m_ii;} + std::vector > m_stoich; - protected: + void addElementaryReaction(const ReactionData& r); + void addThreeBodyReaction(const ReactionData& r); + void addFalloffReaction(const ReactionData& r); - int m_kk, m_nfall; + void installReagents(const ReactionData& r); - array_int m_fallindx; + void installGroups(int irxn, const std::vector& r, + const std::vector& p); + void updateKc(); - Rate1 m_falloff_low_rates; - Rate1 m_falloff_high_rates; - Rate1 m_rates; - - mutable std::map > m_index; - - FalloffMgr m_falloffn; - - ThirdBodyMgr m_3b_concm; - ThirdBodyMgr m_falloff_concm; - - 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; - - GasKineticsData* m_kdata; - - array_fp m_conc; - void processFalloffReactions(); - array_fp m_grt; - - - private: - - int reactionNumber(){ return m_ii;} - std::vector > m_stoich; - - void addElementaryReaction(const ReactionData& r); - void addThreeBodyReaction(const ReactionData& r); - void addFalloffReaction(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; - }; + 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 2ed047f16..451116ada 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -217,10 +217,14 @@ namespace Cantera { * * These routines are basically wrappers around the derived copy * constructor. + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object */ - Kinetics *InterfaceKinetics::duplMyselfAsKinetics() const { - InterfaceKinetics* tp = new InterfaceKinetics(*this); - return dynamic_cast(tp); + Kinetics *InterfaceKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const { + InterfaceKinetics* iK = new InterfaceKinetics(*this); + iK->assignShallowPointers(tpVector); + return dynamic_cast(iK); } //==================================================================================================================== // Update properties that depend on temperature diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 0bb06030c..b6b9208d9 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -115,17 +115,19 @@ namespace Cantera { InterfaceKinetics& operator=(const InterfaceKinetics &right); - //! Duplication routine for objects which inherit from - //! Kinetics + + //! Duplication routine for objects which inherit from Kinetics /*! - * This virtual routine can be used to duplicate %InterfaceKinetics objects + * This virtual routine can be used to duplicate %Kinetics objects * inherited from %Kinetics even if the application only has * a pointer to %Kinetics to work with. * - * These routines are basically wrappers around the derived copy - * constructor. + * These routines are basically wrappers around the derived copy constructor. + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object */ - virtual Kinetics *duplMyselfAsKinetics() const; + virtual Kinetics *duplMyselfAsKinetics(const std::vector & tpVector) const; //! Return the ID of the kinetics object virtual int ID() const; diff --git a/Cantera/src/kinetics/Kinetics.cpp b/Cantera/src/kinetics/Kinetics.cpp index 225ba887d..e935c0410 100644 --- a/Cantera/src/kinetics/Kinetics.cpp +++ b/Cantera/src/kinetics/Kinetics.cpp @@ -32,11 +32,17 @@ namespace Cantera { Kinetics::Kinetics() : m_ii(0), m_nTotalSpecies(0), + m_perturb(0), + m_reactants(0), + m_products(0), m_thermo(0), + m_start(0), + m_phaseindex(), m_index(-1), m_surfphase(-1), m_rxnphase(-1), - m_mindim(4) + m_mindim(4), + m_dummygroups(0) { } @@ -51,16 +57,22 @@ namespace Cantera { Kinetics::Kinetics(const Kinetics &right) : m_ii(0), m_nTotalSpecies(0), + m_perturb(0), + m_reactants(0), + m_products(0), m_thermo(0), + m_start(0), + m_phaseindex(), m_index(-1), m_surfphase(-1), m_rxnphase(-1), - m_mindim(4) + m_mindim(4), + m_dummygroups(0) { /* * Call the assignment operator */ - *this = operator=(right); + *this = right; } // Assignment operator @@ -96,7 +108,7 @@ namespace Cantera { return *this; } - + //==================================================================================================================== // Duplication routine for objects which inherit from // Kinetics /* @@ -107,22 +119,48 @@ namespace Cantera { * These routines are basically wrappers around the derived copy * constructor. */ - Kinetics *Kinetics::duplMyselfAsKinetics() const { - Kinetics* tp = new Kinetics(*this); - return tp; + Kinetics *Kinetics::duplMyselfAsKinetics(const std::vector & tpVector) const { + Kinetics* ko = new Kinetics(*this); + + ko->assignShallowPointers(tpVector); + return ko; } - - - + //==================================================================================================================== int Kinetics::ID() const { return 0; } - + //==================================================================================================================== int Kinetics::type() const { return 0; } + //==================================================================================================================== + void Kinetics::assignShallowPointers(const std::vector & tpVector) { + size_t ns = tpVector.size(); + if (ns != m_thermo.size()) { + throw CanteraError(" Kinetics::assignShallowPointers", + " Number of ThermoPhase objects arent't the same"); + } + for (size_t i = 0; i < ns; i++) { + ThermoPhase *ntp = tpVector[i]; + ThermoPhase *otp = m_thermo[i]; + if (ntp->id() != otp->id()) { + throw CanteraError(" Kinetics::assignShallowPointers", + " id() of the ThermoPhase objects isn't the same"); + } + if (ntp->eosType() != otp->eosType()) { + throw CanteraError(" Kinetics::assignShallowPointers", + " eosType() of the ThermoPhase objects isn't the same"); + } + if (ntp->nSpecies() != otp->nSpecies()) { + throw CanteraError(" Kinetics::assignShallowPointers", + " Number of ThermoPhase objects isn't the same"); + } + m_thermo[i] = tpVector[i]; + } + } + //==================================================================================================================== /** * Takes as input an array of properties for all species in the * mechanism and copies those values beloning to a particular @@ -284,17 +322,13 @@ namespace Cantera { if (type() == cEdgeKinetics) ptype = cEdge; else if (type() == cInterfaceKinetics) ptype = cSurf; if (thermo.eosType() == ptype) { - // if (m_surfphase >= 0) { - // throw CanteraError("Kinetics::addPhase", - // "cannot add more than one surface phase"); - // } m_surfphase = nPhases(); m_rxnphase = nPhases(); } m_thermo.push_back(&thermo); m_phaseindex[m_thermo.back()->id()] = nPhases(); } - + void Kinetics::finalize() { m_nTotalSpecies = 0; int np = nPhases(); @@ -305,9 +339,9 @@ namespace Cantera { } - //! Private function of the class Kinetics, indicating that a function - //! inherited from the base class hasn't had a definition assigned to it - /*! + // Private function of the class Kinetics, indicating that a function + // inherited from the base class hasn't had a definition assigned to it + /* * @param m String message */ void Kinetics::err(std::string m) const { diff --git a/Cantera/src/kinetics/Kinetics.h b/Cantera/src/kinetics/Kinetics.h index f99b0a186..8559d4137 100644 --- a/Cantera/src/kinetics/Kinetics.h +++ b/Cantera/src/kinetics/Kinetics.h @@ -169,18 +169,35 @@ namespace Cantera { Kinetics& operator=(const Kinetics &right); - //! Duplication routine for objects which inherit from - //! Kinetics + //! Duplication routine for objects which inherit from Kinetics /*! * This virtual routine can be used to duplicate %Kinetics objects * inherited from %Kinetics even if the application only has * a pointer to %Kinetics to work with. * - * These routines are basically wrappers around the derived copy - * constructor. + * These routines are basically wrappers around the derived copy constructor. + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object */ - virtual Kinetics *duplMyselfAsKinetics() const; + virtual Kinetics *duplMyselfAsKinetics(const std::vector & tpVector) const; + //! Reassign the shallow pointers within the %FKinetics object + /*! + * This type or routine is absolute necessary because the Kinetics object doesn't + * own the ThermoPhase objects. After a duplication, we need to point to different + * ThermoPhase objects. + * + * We check that the ThermoPhase objects are alligned in the same order and have + * the following identical properties to the ones that they are replacing. + * id() + * eosType() + * nSpecies() + * + * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the + * m_thermo vector within this object + */ + virtual void assignShallowPointers(const std::vector & tpVector); //! Identifies the subclass of the Kinetics manager type. /*! @@ -875,13 +892,13 @@ namespace Cantera { */ //@{ - /// The current value of the multiplier for reaction i. + //! The current value of the multiplier for reaction i. /*! * @param i index of the reaction */ doublereal multiplier(int i) const {return m_perturb[i];} - /// Set the multiplier for reaction i to f. + //! Set the multiplier for reaction i to f. /*! * @param i index of the reaction * @param f value of the multiplier. @@ -965,8 +982,7 @@ namespace Cantera { */ std::vector m_products; - //! m_thermo is a vector of pointers to ThermoPhase - //! objects. + //! m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator /*! * For homogeneous kinetics applications, this vector * will only have one entry. For interfacial reactions, this @@ -999,21 +1015,27 @@ namespace Cantera { * -1. */ std::map m_phaseindex; + //! Index of the Kinetics Manager int m_index; - /** - * Index in the list of phases of the one surface phase. - */ + + //! Index in the list of phases of the one surface phase. + /*! + * + */ int m_surfphase; - /** - * Index in the list of phases of the one phase where the reactions - * occur. + + //! Phase Index where reactions are assumed to be taking place + /*! + * We calculate this by assuming that the phase with the lowest dimensionality is the phase where reactions + * are taking place + * @deprecated */ int m_rxnphase; - /// number of spatial dimensions of lowest-dimensional phase. + //! number of spatial dimensions of lowest-dimensional phase. int m_mindim; private: @@ -1021,9 +1043,11 @@ namespace Cantera { //! Vector of group lists std::vector m_dummygroups; - //! Function for unhandled situations + + //! Private function of the class Kinetics, indicating that a function + //! inherited from the base class hasn't had a definition assigned to it /*! - * @param m String error message + * @param m String message */ void err(std::string m) const;