From 9b90beb0e063b3da8da8ff6c4d7f9ec0f69b5859 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 10 Dec 2009 23:20:44 +0000 Subject: [PATCH] Merged changes from LiquidTransportDevelopment There were two changes: Additions to support some reversible reactions which have fractional stoichiometries. These were not allowed previously. However, if the species with fractional stoichiometries are in single-species phases, this is ok, and actually necessary for some solid-phase reactions. Removed comments from DEBUG_MODE about about too high or low an activation for electrode reactions. --- Cantera/src/kinetics/AqueousKinetics.cpp | 4 +- Cantera/src/kinetics/AqueousKinetics.h | 6 +- Cantera/src/kinetics/EdgeKinetics.h | 6 +- Cantera/src/kinetics/Enhanced3BConc.h | 6 +- Cantera/src/kinetics/FalloffFactory.cpp | 6 +- Cantera/src/kinetics/FalloffFactory.h | 4 +- Cantera/src/kinetics/FalloffMgr.h | 6 +- Cantera/src/kinetics/GasKinetics.h | 6 +- Cantera/src/kinetics/GasKineticsWriter.h | 6 +- Cantera/src/kinetics/Group.cpp | 6 +- Cantera/src/kinetics/Group.h | 6 +- Cantera/src/kinetics/ImplicitChem.cpp | 6 +- Cantera/src/kinetics/ImplicitChem.h | 6 +- Cantera/src/kinetics/ImplicitSurfChem.cpp | 6 +- Cantera/src/kinetics/ImplicitSurfChem.h | 6 +- Cantera/src/kinetics/InterfaceKinetics.cpp | 959 +++++++------ Cantera/src/kinetics/InterfaceKinetics.h | 6 +- Cantera/src/kinetics/Kinetics.cpp | 4 +- Cantera/src/kinetics/Kinetics.h | 4 +- Cantera/src/kinetics/KineticsFactory.cpp | 6 +- Cantera/src/kinetics/KineticsFactory.h | 6 +- Cantera/src/kinetics/Makefile.in | 6 +- Cantera/src/kinetics/RateCoeffMgr.h | 6 +- Cantera/src/kinetics/ReactionData.h | 90 +- Cantera/src/kinetics/ReactionPath.cpp | 6 +- Cantera/src/kinetics/ReactionPath.h | 6 +- Cantera/src/kinetics/ReactionStoichMgr.cpp | 370 +++-- Cantera/src/kinetics/ReactionStoichMgr.h | 6 +- Cantera/src/kinetics/RxnRates.h | 6 +- Cantera/src/kinetics/StoichManager.h | 1509 ++++++++++---------- Cantera/src/kinetics/ThirdBodyMgr.h | 6 +- Cantera/src/kinetics/importKinetics.cpp | 384 ++--- Cantera/src/kinetics/importKinetics.h | 6 +- Cantera/src/kinetics/reaction_defs.h | 6 +- Cantera/src/kinetics/solveSP.cpp | 2 +- Cantera/src/kinetics/solveSP.h | 2 +- 36 files changed, 1782 insertions(+), 1700 deletions(-) diff --git a/Cantera/src/kinetics/AqueousKinetics.cpp b/Cantera/src/kinetics/AqueousKinetics.cpp index 7590bba54..0559dec2b 100644 --- a/Cantera/src/kinetics/AqueousKinetics.cpp +++ b/Cantera/src/kinetics/AqueousKinetics.cpp @@ -11,8 +11,8 @@ * U.S. Government retains certain rights in this software. */ /* - * $Date: 2009/02/11 01:50:54 $ - * $Revision: 1.1 $ + * $Date$ + * $Revision$ */ diff --git a/Cantera/src/kinetics/AqueousKinetics.h b/Cantera/src/kinetics/AqueousKinetics.h index bcdd1d9e7..a23bc0626 100644 --- a/Cantera/src/kinetics/AqueousKinetics.h +++ b/Cantera/src/kinetics/AqueousKinetics.h @@ -3,9 +3,9 @@ * * @ingroup chemkinetics * - * $Author: hkmoffa $ - * $Revision: 1.1 $ - * $Date: 2009/02/11 01:50:57 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/EdgeKinetics.h b/Cantera/src/kinetics/EdgeKinetics.h index c5e2168ca..b19e7b3e6 100644 --- a/Cantera/src/kinetics/EdgeKinetics.h +++ b/Cantera/src/kinetics/EdgeKinetics.h @@ -5,9 +5,9 @@ * @ingroup electrochem */ -/* $Author: hkmoffa $ - * $Revision: 1.4 $ - * $Date: 2008/12/16 20:32:18 $ +/* $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/Enhanced3BConc.h b/Cantera/src/kinetics/Enhanced3BConc.h index 7d927eafd..0aa892a64 100644 --- a/Cantera/src/kinetics/Enhanced3BConc.h +++ b/Cantera/src/kinetics/Enhanced3BConc.h @@ -2,9 +2,9 @@ * @file Enhanced3BConc.h */ -/* $Author: dggoodwin $ - * $Date: 2007/05/04 14:27:23 $ - * $Revision: 1.1 $ +/* $Author$ + * $Date$ + * $Revision$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/FalloffFactory.cpp b/Cantera/src/kinetics/FalloffFactory.cpp index 1e9108df7..fa3f41cf0 100644 --- a/Cantera/src/kinetics/FalloffFactory.cpp +++ b/Cantera/src/kinetics/FalloffFactory.cpp @@ -2,9 +2,9 @@ * @file FalloffFactory.cpp */ -/* $Author: hkmoffa $ - * $Date: 2008/12/29 21:34:08 $ - * $Revision: 1.6 $ +/* $Author$ + * $Date$ + * $Revision$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/FalloffFactory.h b/Cantera/src/kinetics/FalloffFactory.h index fd1605cd7..cbf3a88af 100644 --- a/Cantera/src/kinetics/FalloffFactory.h +++ b/Cantera/src/kinetics/FalloffFactory.h @@ -6,8 +6,8 @@ */ /* - * $Date: 2008/12/29 21:34:08 $ - * $Revision: 1.4 $ + * $Date$ + * $Revision$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/FalloffMgr.h b/Cantera/src/kinetics/FalloffMgr.h index d09cff962..a36d7bcce 100644 --- a/Cantera/src/kinetics/FalloffMgr.h +++ b/Cantera/src/kinetics/FalloffMgr.h @@ -1,9 +1,9 @@ /** * @file FalloffMgr.h * - * $Author: dggoodwin $ - * $Date: 2007/05/10 03:28:32 $ - * $Revision: 1.2 $ + * $Author$ + * $Date$ + * $Revision$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/GasKinetics.h b/Cantera/src/kinetics/GasKinetics.h index bd2c55d0b..7dee80dcd 100644 --- a/Cantera/src/kinetics/GasKinetics.h +++ b/Cantera/src/kinetics/GasKinetics.h @@ -3,9 +3,9 @@ * * @ingroup chemkinetics * - * $Author: hkmoffa $ - * $Revision: 1.8 $ - * $Date: 2008/12/29 21:34:08 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/GasKineticsWriter.h b/Cantera/src/kinetics/GasKineticsWriter.h index fe3be0de0..88b07a213 100644 --- a/Cantera/src/kinetics/GasKineticsWriter.h +++ b/Cantera/src/kinetics/GasKineticsWriter.h @@ -2,9 +2,9 @@ * * @file GasKineticsWriter.h * - * $Author: hkmoffa $ - * $Revision: 1.3 $ - * $Date: 2008/12/17 17:09:37 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/Group.cpp b/Cantera/src/kinetics/Group.cpp index b799def71..4c58727d0 100644 --- a/Cantera/src/kinetics/Group.cpp +++ b/Cantera/src/kinetics/Group.cpp @@ -3,9 +3,9 @@ * * Implementation file for the Group class used in reaction path analysis. * - * $Author: hkmoffa $ - * $Revision: 1.2 $ - * $Date: 2008/12/29 21:34:08 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/Group.h b/Cantera/src/kinetics/Group.h index c1e2ebdca..82e13b74f 100644 --- a/Cantera/src/kinetics/Group.h +++ b/Cantera/src/kinetics/Group.h @@ -1,9 +1,9 @@ /** * @file Group.h * - * $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ + * $Author$ + * $Revision$ + * $Date$ */ diff --git a/Cantera/src/kinetics/ImplicitChem.cpp b/Cantera/src/kinetics/ImplicitChem.cpp index 57cad99db..5a05e2f76 100644 --- a/Cantera/src/kinetics/ImplicitChem.cpp +++ b/Cantera/src/kinetics/ImplicitChem.cpp @@ -2,9 +2,9 @@ * @file ImplicitChem.cpp */ -/* $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ +/* $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ImplicitChem.h b/Cantera/src/kinetics/ImplicitChem.h index 7d17c06b3..96180bdee 100644 --- a/Cantera/src/kinetics/ImplicitChem.h +++ b/Cantera/src/kinetics/ImplicitChem.h @@ -1,9 +1,9 @@ /** * @file ImplicitChem.h * - * $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ImplicitSurfChem.cpp b/Cantera/src/kinetics/ImplicitSurfChem.cpp index 1801b7001..28e82d45a 100644 --- a/Cantera/src/kinetics/ImplicitSurfChem.cpp +++ b/Cantera/src/kinetics/ImplicitSurfChem.cpp @@ -6,9 +6,9 @@ */ /* - * $Author: hkmoffa $ - * $Revision: 1.4 $ - * $Date: 2007/08/29 19:57:48 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ImplicitSurfChem.h b/Cantera/src/kinetics/ImplicitSurfChem.h index 57b233998..adc0fcdce 100644 --- a/Cantera/src/kinetics/ImplicitSurfChem.h +++ b/Cantera/src/kinetics/ImplicitSurfChem.h @@ -6,9 +6,9 @@ */ /* - * $Author: hkmoffa $ - * $Revision: 1.3 $ - * $Date: 2007/08/23 21:43:07 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index 89dcbdd99..98df1c8fd 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -25,45 +25,45 @@ using namespace std; namespace Cantera { - ////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////// - /** - * Construct an empty InterfaceKinetics reaction mechanism. - * @param thermo This is an optional parameter that may be - * used to initialize the inherited Kinetics class with - * one ThermoPhase class object -> in other words it's - * useful for initialization of homogeneous kinetics - * mechanisms. - */ - InterfaceKinetics:: - InterfaceKinetics(thermo_t* thermo) : - Kinetics(), - m_kk(0), - m_redo_rates(false), - m_nirrev(0), - m_nrev(0), - m_surf(0), - m_integrator(0), - m_finalized(false), - m_has_coverage_dependence(false), - m_has_electrochem_rxns(false), - m_ioFlag(0) - { - if (thermo != 0) addPhase(*thermo); - m_kdata = new InterfaceKineticsData; - m_kdata->m_temp = 0.0; - } + /** + * Construct an empty InterfaceKinetics reaction mechanism. + * @param thermo This is an optional parameter that may be + * used to initialize the inherited Kinetics class with + * one ThermoPhase class object -> in other words it's + * useful for initialization of homogeneous kinetics + * mechanisms. + */ + InterfaceKinetics:: + InterfaceKinetics(thermo_t* thermo) : + Kinetics(), + m_kk(0), + m_redo_rates(false), + m_nirrev(0), + m_nrev(0), + m_surf(0), + m_integrator(0), + m_finalized(false), + m_has_coverage_dependence(false), + m_has_electrochem_rxns(false), + m_ioFlag(0) + { + if (thermo != 0) addPhase(*thermo); + m_kdata = new InterfaceKineticsData; + m_kdata->m_temp = 0.0; + } - /** - * Destructor - */ - InterfaceKinetics:: - ~InterfaceKinetics(){ - delete m_kdata; - if (m_integrator) { - delete m_integrator; - } + /** + * Destructor + */ + InterfaceKinetics:: + ~InterfaceKinetics(){ + delete m_kdata; + if (m_integrator) { + delete m_integrator; } + } // Copy Constructor for the %InterfaceKinetics object. /* @@ -154,41 +154,41 @@ namespace Cantera { } - /** - * Update properties that depend on temperature - * - */ - void InterfaceKinetics:: - _update_rates_T() { - _update_rates_phi(); - if (m_has_coverage_dependence) { - m_surf->getCoverages(DATA_PTR(m_conc)); - m_rates.update_C(DATA_PTR(m_conc)); - m_redo_rates = true; - } - doublereal T = thermo(surfacePhaseIndex()).temperature(); - if (T != m_kdata->m_temp || m_redo_rates) { - m_kdata->m_logtemp = log(T); - m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn)); - if (m_has_electrochem_rxns) - applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn)); - m_kdata->m_temp = T; - updateKc(); - m_kdata->m_ROP_ok = false; - m_redo_rates = false; - } + /** + * Update properties that depend on temperature + * + */ + void InterfaceKinetics:: + _update_rates_T() { + _update_rates_phi(); + if (m_has_coverage_dependence) { + m_surf->getCoverages(DATA_PTR(m_conc)); + m_rates.update_C(DATA_PTR(m_conc)); + m_redo_rates = true; } + doublereal T = thermo(surfacePhaseIndex()).temperature(); + if (T != m_kdata->m_temp || m_redo_rates) { + m_kdata->m_logtemp = log(T); + m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn)); + if (m_has_electrochem_rxns) + applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn)); + m_kdata->m_temp = T; + updateKc(); + m_kdata->m_ROP_ok = false; + m_redo_rates = false; + } + } - void InterfaceKinetics:: - _update_rates_phi() { - int np = nPhases(); - for (int n = 0; n < np; n++) { - if (thermo(n).electricPotential() != m_phi[n]) { - m_phi[n] = thermo(n).electricPotential(); - m_redo_rates = true; - } - } + void InterfaceKinetics:: + _update_rates_phi() { + int np = nPhases(); + for (int n = 0; n < np; n++) { + if (thermo(n).electricPotential() != m_phi[n]) { + m_phi[n] = thermo(n).electricPotential(); + m_redo_rates = true; + } } + } /** @@ -265,78 +265,78 @@ namespace Cantera { } - void InterfaceKinetics::checkPartialEquil() { - int i, irxn; - vector_fp dmu(nTotalSpecies(), 0.0); - vector_fp rmu(nReactions(), 0.0); - vector_fp frop(nReactions(), 0.0); - vector_fp rrop(nReactions(), 0.0); - vector_fp netrop(nReactions(), 0.0); - if (m_nrev > 0) { - doublereal rt = GasConstant*thermo(0).temperature(); - cout << "T = " << thermo(0).temperature() << " " << rt << endl; - int n, nsp, k, ik=0; - //doublereal rt = GasConstant*thermo(0).temperature(); - // doublereal rrt = 1.0/rt; - int np = nPhases(); - doublereal delta; - for (n = 0; n < np; n++) { - thermo(n).getChemPotentials(DATA_PTR(dmu) + m_start[n]); - nsp = thermo(n).nSpecies(); - for (k = 0; k < nsp; k++) { - delta = Faraday * m_phi[n] * thermo(n).charge(k); - //cout << thermo(n).speciesName(k) << " " << (delta+dmu[ik])/rt << " " << dmu[ik]/rt << endl; - dmu[ik] += delta; - ik++; - } - } + void InterfaceKinetics::checkPartialEquil() { + int i, irxn; + vector_fp dmu(nTotalSpecies(), 0.0); + vector_fp rmu(nReactions(), 0.0); + vector_fp frop(nReactions(), 0.0); + vector_fp rrop(nReactions(), 0.0); + vector_fp netrop(nReactions(), 0.0); + if (m_nrev > 0) { + doublereal rt = GasConstant*thermo(0).temperature(); + cout << "T = " << thermo(0).temperature() << " " << rt << endl; + int n, nsp, k, ik=0; + //doublereal rt = GasConstant*thermo(0).temperature(); + // doublereal rrt = 1.0/rt; + int np = nPhases(); + doublereal delta; + for (n = 0; n < np; n++) { + thermo(n).getChemPotentials(DATA_PTR(dmu) + m_start[n]); + nsp = thermo(n).nSpecies(); + for (k = 0; k < nsp; k++) { + delta = Faraday * m_phi[n] * thermo(n).charge(k); + //cout << thermo(n).speciesName(k) << " " << (delta+dmu[ik])/rt << " " << dmu[ik]/rt << endl; + dmu[ik] += delta; + ik++; + } + } - // compute Delta mu^ for all reversible reactions - m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(dmu), DATA_PTR(rmu)); - getFwdRatesOfProgress(DATA_PTR(frop)); - getRevRatesOfProgress(DATA_PTR(rrop)); - getNetRatesOfProgress(DATA_PTR(netrop)); - for (i = 0; i < m_nrev; i++) { - irxn = m_revindex[i]; - cout << "Reaction " << reactionString(irxn) - << " " << rmu[irxn]/rt << endl; - printf("%12.6e %12.6e %12.6e %12.6e \n", - frop[irxn], rrop[irxn], netrop[irxn], - netrop[irxn]/(frop[irxn] + rrop[irxn])); - } - } + // compute Delta mu^ for all reversible reactions + m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(dmu), DATA_PTR(rmu)); + getFwdRatesOfProgress(DATA_PTR(frop)); + getRevRatesOfProgress(DATA_PTR(rrop)); + getNetRatesOfProgress(DATA_PTR(netrop)); + for (i = 0; i < m_nrev; i++) { + irxn = m_revindex[i]; + cout << "Reaction " << reactionString(irxn) + << " " << rmu[irxn]/rt << endl; + printf("%12.6e %12.6e %12.6e %12.6e \n", + frop[irxn], rrop[irxn], netrop[irxn], + netrop[irxn]/(frop[irxn] + rrop[irxn])); + } + } + } + + + /** + * Get the equilibrium constants of all reactions, whether + * reversible or not. + */ + void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) { + int i; + + 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++; + } } + fill(kc, kc + m_ii, 0.0); - /** - * Get the equilibrium constants of all reactions, whether - * reversible or not. - */ - void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) { - int i; + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc); - 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++; - } - } - - fill(kc, kc + m_ii, 0.0); - - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc); - - for (i = 0; i < m_ii; i++) { - kc[i] = exp(-kc[i]*rrt); - } + for (i = 0; i < m_ii; i++) { + kc[i] = exp(-kc[i]*rrt); } + } // Returns the Species creation rates [kmol/m^2/s]. /* @@ -414,10 +414,17 @@ namespace Cantera { DATA_PTR(m_rwork)); // Modify the reaction rates. Only modify those with a - // non-zero activation energy, and do not decrease the - // activation energy below zero. + // non-zero activation energy. Below we decrease the + // activation energy below zero but in some debug modes + // we print out a warning message about this. + /* + * NOTE, there is some discussion about this point. + * Should we decrease the activiation energy below zero? + * I don't think this has been decided in any definative way. + * The treatment below is numerically more stable, however. + */ doublereal eamod; -#ifdef DEBUG_MODE +#ifdef DEBUG_KIN_MODE double ea; #endif int nct = m_beta.size(); @@ -426,7 +433,7 @@ namespace Cantera { irxn = m_ctrxn[i]; eamod = m_beta[i]*m_rwork[irxn]; if (eamod != 0.0 && m_E[irxn] != 0.0) { -#ifdef DEBUG_MODE +#ifdef DEBUG_KIN_MODE ea = GasConstant * m_E[irxn]; if (eamod + ea < 0.0) { writelog("Warning: act energy mod too large!\n"); @@ -447,162 +454,162 @@ namespace Cantera { - /** - * Update the rates of progress of the reactions in the reaciton - * mechanism. This routine operates on internal data. - */ - void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) { + /** + * Update the rates of progress of the reactions in the reaciton + * mechanism. This routine operates on internal data. + */ + void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) { - // _update_rates_T(); - // _update_rates_C(); - updateROP(); + // _update_rates_T(); + // _update_rates_C(); + updateROP(); - const vector_fp& rf = m_kdata->m_rfn; + const vector_fp& rf = m_kdata->m_rfn; - // copy rate coefficients into kfwd - copy(rf.begin(), rf.end(), kfwd); + // copy rate coefficients into kfwd + copy(rf.begin(), rf.end(), kfwd); - // multiply by perturbation factor - multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin()); + // multiply by perturbation factor + multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin()); + } + + + /** + * Update the rates of progress of the reactions in the reaciton + * mechanism. This routine operates on internal data. + */ + void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) { + getFwdRateConstants(krev); + if (doIrreversible) { + doublereal *tmpKc = DATA_PTR(m_kdata->m_ropnet); + getEquilibriumConstants(tmpKc); + for (int i = 0; i < m_ii; i++) { + krev[i] /= tmpKc[i]; + } } - - - /** - * Update the rates of progress of the reactions in the reaciton - * mechanism. This routine operates on internal data. - */ - void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) { - getFwdRateConstants(krev); - if (doIrreversible) { - doublereal *tmpKc = DATA_PTR(m_kdata->m_ropnet); - getEquilibriumConstants(tmpKc); - for (int i = 0; i < m_ii; i++) { - krev[i] /= tmpKc[i]; - } - } - else { - const vector_fp& rkc = m_kdata->m_rkcn; - multiply_each(krev, krev + nReactions(), rkc.begin()); - } + else { + const vector_fp& rkc = m_kdata->m_rkcn; + multiply_each(krev, krev + nReactions(), rkc.begin()); } + } - void InterfaceKinetics::getActivationEnergies(doublereal *E) { - copy(m_E.begin(), m_E.end(), E); - } + void InterfaceKinetics::getActivationEnergies(doublereal *E) { + copy(m_E.begin(), m_E.end(), E); + } - /** - * Update the rates of progress of the reactions in the reaction - * mechanism. This routine operates on internal data. - */ - void InterfaceKinetics::updateROP() { + /** + * Update the rates of progress of the reactions in the reaction + * mechanism. This routine operates on internal data. + */ + void InterfaceKinetics::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 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(DATA_PTR(m_conc), DATA_PTR(ropf)); - //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); + // multiply ropf by concentration products + m_rxnstoich.multiplyReactants(DATA_PTR(m_conc), DATA_PTR(ropf)); + //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); - // for reversible reactions, multiply ropr by concentration - // products - m_rxnstoich.multiplyRevProducts(DATA_PTR(m_conc), - DATA_PTR(ropr)); - //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); + // for reversible reactions, multiply ropr by concentration + // products + m_rxnstoich.multiplyRevProducts(DATA_PTR(m_conc), + DATA_PTR(ropr)); + //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); - // do global reactions - //m_globalReactantStoich.power(m_conc.begin(), ropf.begin()); + // do global reactions + //m_globalReactantStoich.power(m_conc.begin(), ropf.begin()); - for (int j = 0; j != m_ii; ++j) { - ropnet[j] = ropf[j] - ropr[j]; - } + for (int j = 0; j != m_ii; ++j) { + ropnet[j] = ropf[j] - ropr[j]; + } - m_kdata->m_ROP_ok = true; - } + m_kdata->m_ROP_ok = true; + } - /** - * - * 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 + /** + * + * 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 InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) { + /* + * Get the chemical potentials of the species in the + * ideal gas solution. */ - void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) { - /* - * Get the chemical potentials of the species in the - * ideal gas solution. - */ - int np = nPhases(); - int n; - for (n = 0; n < np; n++) { - thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]); - } - //for (n = 0; n < m_grt.size(); n++) { - // cout << n << "G_RT = " << m_grt[n] << endl; - //} - - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); + int np = nPhases(); + int n; + for (n = 0; n < np; n++) { + thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]); } + //for (n = 0; n < m_grt.size(); n++) { + // cout << n << "G_RT = " << m_grt[n] << endl; + //} + + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), 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 + /** + * + * 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 InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) { + /* + * Get the partial molar enthalpy of all species in the + * ideal gas. */ - void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) { - /* - * Get the partial molar enthalpy of all species in the - * ideal gas. - */ - int np = nPhases(); - int n; - for (n = 0; n < np; n++) { - thermo(n).getPartialMolarEnthalpies(DATA_PTR(m_grt) + m_start[n]); - } - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); + int np = nPhases(); + int n; + for (n = 0; n < np; n++) { + thermo(n).getPartialMolarEnthalpies(DATA_PTR(m_grt) + m_start[n]); } + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); + } // Return the vector of values for the change in @@ -634,160 +641,160 @@ namespace Cantera { m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), 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 + /** + * + * 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 InterfaceKinetics::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. */ - void InterfaceKinetics::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. - */ - int np = nPhases(); - int n; - for (n = 0; n < np; n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_grt) + m_start[n]); - } - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); + int np = nPhases(); + int n; + for (n = 0; n < np; n++) { + thermo(n).getStandardChemPotentials(DATA_PTR(m_grt) + m_start[n]); } - - /** - * - * 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 + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. */ - void InterfaceKinetics::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. - */ - int np = nPhases(); - int n; - for (n = 0; n < np; n++) { - thermo(n).getEnthalpy_RT(DATA_PTR(m_grt) + m_start[n]); - } - 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, DATA_PTR(m_grt), deltaH); - } + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); + } - /********************************************************************* - * - * 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 + /** + * + * 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 InterfaceKinetics::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. */ - void InterfaceKinetics::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. - */ - int np = nPhases(); - int n; - for (n = 0; n < np; n++) { - thermo(n).getEntropy_R(DATA_PTR(m_grt) + m_start[n]); - } - 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, DATA_PTR(m_grt), deltaS); + int np = nPhases(); + int n; + for (n = 0; n < np; n++) { + thermo(n).getEnthalpy_RT(DATA_PTR(m_grt) + m_start[n]); } - - - /** - * Add a single reaction to the mechanism. This routine - * must be called after init() and before finalize(). - * This function branches on the types of reactions allowed - * by the interfaceKinetics manager in order to install - * the reaction correctly in the manager. - * The manager allows the following reaction types - * Elementary - * Surface - * Global - * There is no difference between elementary and surface - * reactions. + 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. */ - void InterfaceKinetics:: - addReaction(const ReactionData& r) { + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH); + } - addElementaryReaction(r); + /********************************************************************* + * + * 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 InterfaceKinetics::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. + */ + int np = nPhases(); + int n; + for (n = 0; n < np; n++) { + thermo(n).getEntropy_R(DATA_PTR(m_grt) + m_start[n]); + } + 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, DATA_PTR(m_grt), deltaS); + } - // operations common to all reaction types - installReagents( r ); - //installGroups(reactionNumber(), r.rgroups, r.pgroups); - incrementRxnCount(); - m_rxneqn.push_back(r.equation); + + /** + * Add a single reaction to the mechanism. This routine + * must be called after init() and before finalize(). + * This function branches on the types of reactions allowed + * by the interfaceKinetics manager in order to install + * the reaction correctly in the manager. + * The manager allows the following reaction types + * Elementary + * Surface + * Global + * There is no difference between elementary and surface + * reactions. + */ + void InterfaceKinetics:: + addReaction(const ReactionData& r) { + + addElementaryReaction(r); + + // operations common to all reaction types + installReagents( r ); + //installGroups(reactionNumber(), r.rgroups, r.pgroups); + incrementRxnCount(); + m_rxneqn.push_back(r.equation); + } + + + void InterfaceKinetics:: + addElementaryReaction(const ReactionData& r) { + int iloc; + + // install rate coeff calculator + + vector_fp rp = r.rateCoeffParameters; + int ncov = r.cov.size(); + if (ncov > 3) { + m_has_coverage_dependence = true; + } + for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); + + iloc = m_rates.install( reactionNumber(), + r.rateCoeffType, rp.size(), + DATA_PTR(rp) ); + + // store activation energy + m_E.push_back(r.rateCoeffParameters[2]); + + if (r.beta > 0.0) { + m_has_electrochem_rxns = true; + m_beta.push_back(r.beta); + m_ctrxn.push_back(reactionNumber()); } - - void InterfaceKinetics:: - addElementaryReaction(const ReactionData& r) { - int iloc; - - // install rate coeff calculator - - vector_fp rp = r.rateCoeffParameters; - int ncov = r.cov.size(); - if (ncov > 3) { - m_has_coverage_dependence = true; - } - for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); - - iloc = m_rates.install( reactionNumber(), - r.rateCoeffType, rp.size(), - DATA_PTR(rp) ); - - // store activation energy - m_E.push_back(r.rateCoeffParameters[2]); - - if (r.beta > 0.0) { - m_has_electrochem_rxns = true; - m_beta.push_back(r.beta); - m_ctrxn.push_back(reactionNumber()); - } - - // add constant term to rate coeff value vector - m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); - registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc); - } + // add constant term to rate coeff value vector + m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); + registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc); + } void InterfaceKinetics::setIOFlag(int ioFlag) { @@ -797,32 +804,32 @@ namespace Cantera { } } -// void InterfaceKinetics:: -// addGlobalReaction(const ReactionData& r) { + // void InterfaceKinetics:: + // addGlobalReaction(const ReactionData& r) { -// int iloc; -// // install rate coeff calculator -// vector_fp rp = r.rateCoeffParameters; -// int ncov = r.cov.size(); -// for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); -// iloc = m_rates.install( reactionNumber(), -// r.rateCoeffType, rp.size(), -// rp.begin() ); -// // store activation energy -// m_E.push_back(r.rateCoeffParameters[2]); -// // add constant term to rate coeff value vector -// m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); + // int iloc; + // // install rate coeff calculator + // vector_fp rp = r.rateCoeffParameters; + // int ncov = r.cov.size(); + // for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); + // iloc = m_rates.install( reactionNumber(), + // r.rateCoeffType, rp.size(), + // rp.begin() ); + // // store activation energy + // m_E.push_back(r.rateCoeffParameters[2]); + // // add constant term to rate coeff value vector + // m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); -// int nr = r.order.size(); -// vector_fp ordr(nr); -// for (int n = 0; n < nr; n++) { -// ordr[n] = r.order[n] - r.rstoich[n]; -// } -// m_globalReactantStoich.add( reactionNumber(), -// r.reactants, ordr); + // int nr = r.order.size(); + // vector_fp ordr(nr); + // for (int n = 0; n < nr; n++) { + // ordr[n] = r.order[n] - r.rstoich[n]; + // } + // m_globalReactantStoich.add( reactionNumber(), + // r.reactants, ordr); -// registerReaction( reactionNumber(), GLOBAL_RXN, iloc); -// } + // registerReaction( reactionNumber(), GLOBAL_RXN, iloc); + // } void InterfaceKinetics::installReagents(const ReactionData& r) { @@ -918,50 +925,50 @@ namespace Cantera { } } - /** - * Prepare the class for the addition of reactions. This function - * must be called after instantiation of the class, but before - * any reactions are actually added to the mechanism. - * This function calculates m_kk the number of species in all - * phases participating in the reaction mechanism. We don't know - * m_kk previously, before all phases have been added. - */ - void InterfaceKinetics::init() { - int n; - m_kk = 0; - int np = nPhases(); - for (n = 0; n < np; n++) { - m_kk += thermo(n).nSpecies(); - } - m_rrxn.resize(m_kk); - m_prxn.resize(m_kk); - m_conc.resize(m_kk); - m_mu0.resize(m_kk); - m_grt.resize(m_kk); - m_pot.resize(m_kk, 0.0); - m_phi.resize(np, 0.0); + /** + * Prepare the class for the addition of reactions. This function + * must be called after instantiation of the class, but before + * any reactions are actually added to the mechanism. + * This function calculates m_kk the number of species in all + * phases participating in the reaction mechanism. We don't know + * m_kk previously, before all phases have been added. + */ + void InterfaceKinetics::init() { + int n; + m_kk = 0; + int np = nPhases(); + for (n = 0; n < np; n++) { + m_kk += thermo(n).nSpecies(); } + m_rrxn.resize(m_kk); + m_prxn.resize(m_kk); + m_conc.resize(m_kk); + m_mu0.resize(m_kk); + m_grt.resize(m_kk); + m_pot.resize(m_kk, 0.0); + m_phi.resize(np, 0.0); + } - /** - * Finish adding reactions and prepare for use. This function - * must be called after all reactions are entered into the mechanism - * and before the mechanism is used to calculate reaction rates. - * - * Here, we resize work arrays based on the number of reactions, - * since we don't know this number up to now. - */ - void InterfaceKinetics::finalize() { - m_rwork.resize(nReactions()); - int ks = reactionPhaseIndex(); - if (ks < 0) throw CanteraError("InterfaceKinetics::finalize", - "no surface phase is present."); - m_surf = (SurfPhase*)&thermo(ks); - if (m_surf->nDim() != 2) - throw CanteraError("InterfaceKinetics::finalize", - "expected interface dimension = 2, but got dimension = " - +int2str(m_surf->nDim())); - m_finalized = true; - } + /** + * Finish adding reactions and prepare for use. This function + * must be called after all reactions are entered into the mechanism + * and before the mechanism is used to calculate reaction rates. + * + * Here, we resize work arrays based on the number of reactions, + * since we don't know this number up to now. + */ + void InterfaceKinetics::finalize() { + m_rwork.resize(nReactions()); + int ks = reactionPhaseIndex(); + if (ks < 0) throw CanteraError("InterfaceKinetics::finalize", + "no surface phase is present."); + m_surf = (SurfPhase*)&thermo(ks); + if (m_surf->nDim() != 2) + throw CanteraError("InterfaceKinetics::finalize", + "expected interface dimension = 2, but got dimension = " + +int2str(m_surf->nDim())); + m_finalized = true; + } doublereal InterfaceKinetics::electrochem_beta(int irxn) const{ @@ -975,9 +982,9 @@ namespace Cantera { } - bool InterfaceKinetics::ready() const { - return (m_finalized); - } + bool InterfaceKinetics::ready() const { + return (m_finalized); + } // Advance the surface coverages in time /* diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 261947bf9..3de46a1f1 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -4,9 +4,9 @@ * @ingroup chemkinetics */ /* - * $Author: hkmoffa $ - * $Revision: 1.11 $ - * $Date: 2008/12/17 17:09:37 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/Kinetics.cpp b/Cantera/src/kinetics/Kinetics.cpp index f7ac93ff7..38ea75fe0 100644 --- a/Cantera/src/kinetics/Kinetics.cpp +++ b/Cantera/src/kinetics/Kinetics.cpp @@ -7,8 +7,8 @@ * Kinetics managers calculate rates of progress of species due to homogeneous or heterogeneous kinetics. */ /* - * $Date: 2008/12/16 20:32:18 $ - * $Revision: 1.3 $ + * $Date$ + * $Revision$ */ // Copyright 2001-2004 California Institute of Technology diff --git a/Cantera/src/kinetics/Kinetics.h b/Cantera/src/kinetics/Kinetics.h index 69c6ebca8..25ca5305f 100644 --- a/Cantera/src/kinetics/Kinetics.h +++ b/Cantera/src/kinetics/Kinetics.h @@ -4,8 +4,8 @@ * module documentation (see \ref kineticsmgr and class * \link Cantera::Kinetics Kinetics\endlink). * - * $Date: 2008/12/16 20:32:18 $ - * $Revision: 1.8 $ + * $Date$ + * $Revision$ */ // Copyright 2001-2004 California Institute of Technology diff --git a/Cantera/src/kinetics/KineticsFactory.cpp b/Cantera/src/kinetics/KineticsFactory.cpp index 142bdc36f..0d811f63e 100644 --- a/Cantera/src/kinetics/KineticsFactory.cpp +++ b/Cantera/src/kinetics/KineticsFactory.cpp @@ -3,9 +3,9 @@ */ /* - * $Author: hkmoffa $ - * $Revision: 1.4 $ - * $Date: 2009/02/11 20:01:45 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/KineticsFactory.h b/Cantera/src/kinetics/KineticsFactory.h index 4fdc9c6be..f4336c96d 100644 --- a/Cantera/src/kinetics/KineticsFactory.h +++ b/Cantera/src/kinetics/KineticsFactory.h @@ -3,9 +3,9 @@ */ /* - * $Author: hkmoffa $ - * $Revision: 1.3 $ - * $Date: 2009/02/11 01:50:58 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/Makefile.in b/Cantera/src/kinetics/Makefile.in index 19700d7d6..cfd625c43 100644 --- a/Cantera/src/kinetics/Makefile.in +++ b/Cantera/src/kinetics/Makefile.in @@ -1,8 +1,8 @@ #/bin/sh ############################################################### -# $Author: hkmoffa $ -# $Date: 2009/02/11 01:50:58 $ -# $Revision: 1.10 $ +# $Author$ +# $Date$ +# $Revision$ # # Copyright 2002 California Institute of Technology # diff --git a/Cantera/src/kinetics/RateCoeffMgr.h b/Cantera/src/kinetics/RateCoeffMgr.h index fba1ae9d4..792c14617 100644 --- a/Cantera/src/kinetics/RateCoeffMgr.h +++ b/Cantera/src/kinetics/RateCoeffMgr.h @@ -3,9 +3,9 @@ */ /* - * $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ReactionData.h b/Cantera/src/kinetics/ReactionData.h index 6c7665a31..117f650ef 100644 --- a/Cantera/src/kinetics/ReactionData.h +++ b/Cantera/src/kinetics/ReactionData.h @@ -3,9 +3,9 @@ * */ /* - * $Author: hkmoffa $ - * $Revision: 1.2 $ - * $Date: 2007/06/04 23:05:07 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology @@ -18,48 +18,52 @@ namespace Cantera { - class ReactionData { - public: - ReactionData() { - reactionType = ELEMENTARY_RXN; - number = 0; - rxn_number = 0; - reversible = true; - rateCoeffType = ARRHENIUS; - falloffType = NONE; - error = 0; - equation = ""; - default_3b_eff = 1.0; - global = false; - beta = 0.0; - } - virtual ~ReactionData(){} + class ReactionData { + public: + ReactionData() { + reactionType = ELEMENTARY_RXN; + number = 0; + rxn_number = 0; + reversible = true; + rateCoeffType = ARRHENIUS; + falloffType = NONE; + error = 0; + equation = ""; + default_3b_eff = 1.0; + global = false; + isReversibleWithFrac = false; + beta = 0.0; + } + ~ReactionData(){} - int reactionType; - int number, rxn_number; - vector_int reactants; - vector_int products; - vector_fp order; - vector_fp rstoich; - vector_fp pstoich; - std::vector rgroups; - std::vector pgroups; - std::map thirdBodyEfficiencies; + int reactionType; + int number; + int rxn_number; + vector_int reactants; + vector_int products; + vector_fp rorder; + vector_fp porder; + vector_fp rstoich; + vector_fp pstoich; + std::vector rgroups; + std::vector pgroups; + std::map thirdBodyEfficiencies; - //! True if the current reaction is reversible. False otherwise - bool reversible; - int rateCoeffType; - vector_fp rateCoeffParameters; - vector_fp auxRateCoeffParameters; - int falloffType; - vector_fp falloffParameters; - int error; - std::string equation; - doublereal default_3b_eff; - vector_fp cov; - bool global; - doublereal beta; // for electrochemical reactions - }; + //! True if the current reaction is reversible. False otherwise + bool reversible; + int rateCoeffType; + vector_fp rateCoeffParameters; + vector_fp auxRateCoeffParameters; + int falloffType; + vector_fp falloffParameters; + int error; + std::string equation; + doublereal default_3b_eff; + vector_fp cov; + bool global; + bool isReversibleWithFrac; + doublereal beta; // for electrochemical reactions + }; } #endif diff --git a/Cantera/src/kinetics/ReactionPath.cpp b/Cantera/src/kinetics/ReactionPath.cpp index 4b2feb41f..af42f38a2 100644 --- a/Cantera/src/kinetics/ReactionPath.cpp +++ b/Cantera/src/kinetics/ReactionPath.cpp @@ -4,9 +4,9 @@ */ /* - * $Author: dggoodwin $ - * $Revision: 1.2 $ - * $Date: 2008/02/05 23:36:12 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ReactionPath.h b/Cantera/src/kinetics/ReactionPath.h index 5100648b7..f2ccf67d6 100644 --- a/Cantera/src/kinetics/ReactionPath.h +++ b/Cantera/src/kinetics/ReactionPath.h @@ -3,9 +3,9 @@ * * Classes for reaction path analysis. * - * $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/ReactionStoichMgr.cpp b/Cantera/src/kinetics/ReactionStoichMgr.cpp index dca9fc3b5..6009a1164 100644 --- a/Cantera/src/kinetics/ReactionStoichMgr.cpp +++ b/Cantera/src/kinetics/ReactionStoichMgr.cpp @@ -5,9 +5,9 @@ /// //------------------------------------------------ -// $Author: dggoodwin $ -// $Revision: 1.1 $ -// $Date: 2007/05/04 15:48:44 $ +// $Author$ +// $Revision$ +// $Date$ // turn off warnings under Windows #ifdef WIN32 @@ -67,216 +67,214 @@ namespace Cantera { void ReactionStoichMgr:: add(int rxn, const ReactionData& r) { - vector_int rk; - doublereal frac; - bool isfrac = false; - int n, ns, m, nr = r.reactants.size(); - for (n = 0; n < nr; n++) { - ns = int(r.rstoich[n]); - frac = r.rstoich[n] - 1.0*int(r.rstoich[n]); - if (frac != 0.0) isfrac = true; - for (m = 0; m < ns; m++) { - rk.push_back(r.reactants[n]); - } - } + vector_int rk; + doublereal frac; + bool isfrac = false; + int n, ns, m, nr = r.reactants.size(); + for (n = 0; n < nr; n++) { + ns = int(r.rstoich[n]); + frac = r.rstoich[n] - 1.0*int(r.rstoich[n]); + if (frac != 0.0) isfrac = true; + for (m = 0; m < ns; m++) { + rk.push_back(r.reactants[n]); + } + } - // if the reaction has fractional stoichiometric coefficients - // or specified reaction orders, then add it in a ma - if (isfrac || r.global || rk.size() > 3) { - m_reactants->add(rxn, r.reactants, r.order, r.rstoich); + // if the reaction has fractional stoichiometric coefficients + // or specified reaction orders, then add it in a general reaction + if (isfrac || r.global || rk.size() > 3) { + m_reactants->add(rxn, r.reactants, r.rorder, r.rstoich); #ifdef INCL_STOICH_WRITER - if (m_rwriter) m_rwriter->add(rxn, r.reactants, r.order, r.rstoich); + if (m_rwriter) m_rwriter->add(rxn, r.reactants, r.order, r.rstoich); #endif - } - else { - m_reactants->add( rxn, rk); + } + else { + m_reactants->add( rxn, rk); #ifdef INCL_STOICH_WRITER - if (m_rwriter) m_rwriter->add(rxn, rk); + if (m_rwriter) m_rwriter->add(rxn, rk); #endif - } + } - vector_int pk; - isfrac = false; - int np = r.products.size(); - for (n = 0; n < np; n++) { - ns = int(r.pstoich[n]); - frac = r.pstoich[n] - 1.0*int(r.pstoich[n]); - if (frac != 0.0) isfrac = true; - for (m = 0; m < ns; m++) { - pk.push_back(r.products[n]); - } - } + vector_int pk; + isfrac = false; + int np = r.products.size(); + for (n = 0; n < np; n++) { + ns = int(r.pstoich[n]); + frac = r.pstoich[n] - 1.0*int(r.pstoich[n]); + if (frac != 0.0) isfrac = true; + for (m = 0; m < ns; m++) { + pk.push_back(r.products[n]); + } + } - if (r.reversible) { - if (isfrac) { - throw CanteraError("ReactionStoichMgr::add", - "fractional product stoichiometric coefficients only allowed " - "\nfor irreversible reactions"); - } - if (pk.size() > 3) { - // mod dgg 10/31/06 - m_revproducts->add(rxn, r.products, r.pstoich, r.pstoich); - //m_revproducts->add(rxn, r.products, m_dummy, r.pstoich); - } - else { - m_revproducts->add(rxn, pk); - } - } - else if (isfrac || pk.size() > 3) { - m_irrevproducts->add(rxn, r.products, m_dummy, r.pstoich); - } - else { - m_irrevproducts->add(rxn, pk); - } + if (r.reversible) { + if (isfrac && !r.isReversibleWithFrac) { + throw CanteraError("ReactionStoichMgr::add", + "Fractional product stoichiometric coefficients only allowed " + "\nfor irreversible reactions and most reversible reactions"); + } + if (pk.size() > 3 || r.isReversibleWithFrac) { + m_revproducts->add(rxn, r.products, r.porder, r.pstoich); + } + else { + m_revproducts->add(rxn, pk); + } + } + else if (isfrac || pk.size() > 3) { + m_irrevproducts->add(rxn, r.products, r.porder, r.pstoich); + } + else { + m_irrevproducts->add(rxn, pk); + } } void ReactionStoichMgr:: getCreationRates(int nsp, const doublereal* ropf, - const doublereal* ropr, doublereal* c) { - // zero out the output array - fill(c, c + nsp, 0.0); + const doublereal* ropr, doublereal* c) { + // zero out the output array + fill(c, c + nsp, 0.0); - // the forward direction creates product species - m_revproducts->incrementSpecies(ropf, c); - m_irrevproducts->incrementSpecies(ropf, c); + // the forward direction creates product species + m_revproducts->incrementSpecies(ropf, c); + m_irrevproducts->incrementSpecies(ropf, c); - // the reverse direction creates reactant species - m_reactants->incrementSpecies(ropr, c); + // the reverse direction creates reactant species + m_reactants->incrementSpecies(ropr, c); } - void ReactionStoichMgr:: - getDestructionRates(int nsp, const doublereal* ropf, - const doublereal* ropr, doublereal* d) { - fill(d, d + nsp, 0.0); - // the reverse direction destroys products in reversible reactions - m_revproducts->incrementSpecies(ropr, d); - // the forward direction destroys reactants - m_reactants->incrementSpecies(ropf, d); - } + void ReactionStoichMgr:: + getDestructionRates(int nsp, const doublereal* ropf, + const doublereal* ropr, doublereal* d) { + fill(d, d + nsp, 0.0); + // the reverse direction destroys products in reversible reactions + m_revproducts->incrementSpecies(ropr, d); + // the forward direction destroys reactants + m_reactants->incrementSpecies(ropf, d); + } - void ReactionStoichMgr:: - getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w) { - fill(w, w + nsp, 0.0); - // products are created for positive net rate of progress - m_revproducts->incrementSpecies(ropnet, w); - m_irrevproducts->incrementSpecies(ropnet, w); - // reactants are destroyed for positive net rate of progress - m_reactants->decrementSpecies(ropnet, w); - } + void ReactionStoichMgr:: + getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w) { + fill(w, w + nsp, 0.0); + // products are created for positive net rate of progress + m_revproducts->incrementSpecies(ropnet, w); + m_irrevproducts->incrementSpecies(ropnet, w); + // reactants are destroyed for positive net rate of progress + m_reactants->decrementSpecies(ropnet, w); + } - void ReactionStoichMgr:: - getReactionDelta(int nr, const doublereal* g, doublereal* dg) { - fill(dg, dg + nr, 0.0); - // products add - m_revproducts->incrementReactions(g, dg); - m_irrevproducts->incrementReactions(g, dg); - // reactants subtract - m_reactants->decrementReactions(g, dg); - } + void ReactionStoichMgr:: + getReactionDelta(int nr, const doublereal* g, doublereal* dg) { + fill(dg, dg + nr, 0.0); + // products add + m_revproducts->incrementReactions(g, dg); + m_irrevproducts->incrementReactions(g, dg); + // reactants subtract + m_reactants->decrementReactions(g, dg); + } - void ReactionStoichMgr:: - getRevReactionDelta(int nr, const doublereal* g, doublereal* dg) { - fill(dg, dg + nr, 0.0); - m_revproducts->incrementReactions(g, dg); - m_reactants->decrementReactions(g, dg); - } + void ReactionStoichMgr:: + getRevReactionDelta(int nr, const doublereal* g, doublereal* dg) { + fill(dg, dg + nr, 0.0); + m_revproducts->incrementReactions(g, dg); + m_reactants->decrementReactions(g, dg); + } - void ReactionStoichMgr:: - multiplyReactants(const doublereal* c, doublereal* r) { - m_reactants->multiply(c, r); - } + void ReactionStoichMgr:: + multiplyReactants(const doublereal* c, doublereal* r) { + m_reactants->multiply(c, r); + } - void ReactionStoichMgr:: - multiplyRevProducts(const doublereal* c, doublereal* r) { - m_revproducts->multiply(c, r); - } + void ReactionStoichMgr:: + multiplyRevProducts(const doublereal* c, doublereal* r) { + m_revproducts->multiply(c, r); + } - void ReactionStoichMgr:: - write(string filename) { - ofstream f(filename.c_str()); - f << "namespace mech {" << endl; - writeCreationRates(f); - writeDestructionRates(f); - writeNetProductionRates(f); - writeMultiplyReactants(f); - writeMultiplyRevProducts(f); - f << "} // namespace mech" << endl; - f.close(); + void ReactionStoichMgr:: + write(string filename) { + ofstream f(filename.c_str()); + f << "namespace mech {" << endl; + writeCreationRates(f); + writeDestructionRates(f); + writeNetProductionRates(f); + writeMultiplyReactants(f); + writeMultiplyRevProducts(f); + f << "} // namespace mech" << endl; + f.close(); + } + + void ReactionStoichMgr:: + writeCreationRates(ostream& f) { + f << " void getCreationRates(const doublereal* rf, const doublereal* rb," << endl; + f << " doublereal* c) {" << endl; + map out; + m_revproducts->writeIncrementSpecies("rf",out); + m_irrevproducts->writeIncrementSpecies("rf",out); + m_reactants->writeIncrementSpecies("rb",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = wrapString(b->second); + rhs[1] = '='; + f << " c[" << b->first << "] " << rhs << ";" << endl; } + f << " }" << endl << endl << endl; + } - void ReactionStoichMgr:: - writeCreationRates(ostream& f) { - f << " void getCreationRates(const doublereal* rf, const doublereal* rb," << endl; - f << " doublereal* c) {" << endl; - map out; - m_revproducts->writeIncrementSpecies("rf",out); - m_irrevproducts->writeIncrementSpecies("rf",out); - m_reactants->writeIncrementSpecies("rb",out); - map::iterator b; - for (b = out.begin(); b != out.end(); ++b) { - string rhs = wrapString(b->second); - rhs[1] = '='; - f << " c[" << b->first << "] " << rhs << ";" << endl; - } - f << " }" << endl << endl << endl; - } - - void ReactionStoichMgr:: - writeDestructionRates(ostream& f) { - f << " void getDestructionRates(const doublereal* rf, const doublereal* rb," << endl; - f << " doublereal* d) {" << endl; - map out; - m_revproducts->writeIncrementSpecies("rb",out); - m_reactants->writeIncrementSpecies("rf",out); - map::iterator b; - for (b = out.begin(); b != out.end(); ++b) { - string rhs = wrapString(b->second); - rhs[1] = '='; - f << " d[" << b->first << "] " << rhs << ";" << endl; - } - f << " }" << endl << endl << endl; + void ReactionStoichMgr:: + writeDestructionRates(ostream& f) { + f << " void getDestructionRates(const doublereal* rf, const doublereal* rb," << endl; + f << " doublereal* d) {" << endl; + map out; + m_revproducts->writeIncrementSpecies("rb",out); + m_reactants->writeIncrementSpecies("rf",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = wrapString(b->second); + rhs[1] = '='; + f << " d[" << b->first << "] " << rhs << ";" << endl; } + f << " }" << endl << endl << endl; + } - void ReactionStoichMgr:: - writeNetProductionRates(ostream& f) { - f << " void getNetProductionRates(const doublereal* r, doublereal* w) {" << endl; - map out; - m_revproducts->writeIncrementSpecies("r",out); - m_irrevproducts->writeIncrementSpecies("r",out); - m_reactants->writeDecrementSpecies("r",out); - map::iterator b; - for (b = out.begin(); b != out.end(); ++b) { - string rhs = wrapString(b->second); - rhs[1] = '='; - f << " w[" << b->first << "] " << rhs << ";" << endl; - } - f << " }" << endl << endl << endl; + void ReactionStoichMgr:: + writeNetProductionRates(ostream& f) { + f << " void getNetProductionRates(const doublereal* r, doublereal* w) {" << endl; + map out; + m_revproducts->writeIncrementSpecies("r",out); + m_irrevproducts->writeIncrementSpecies("r",out); + m_reactants->writeDecrementSpecies("r",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = wrapString(b->second); + rhs[1] = '='; + f << " w[" << b->first << "] " << rhs << ";" << endl; } + f << " }" << endl << endl << endl; + } - void ReactionStoichMgr:: - writeMultiplyReactants(ostream& f) { - f << " void multiplyReactants(const doublereal* c, doublereal* r) {" << endl; - map out; - m_reactants->writeMultiply("c",out); - map::iterator b; - for (b = out.begin(); b != out.end(); ++b) { - string rhs = b->second; - f << " r[" << b->first << "] *= " << rhs << ";" << endl; - } - f << " }" << endl << endl << endl; - } + void ReactionStoichMgr:: + writeMultiplyReactants(ostream& f) { + f << " void multiplyReactants(const doublereal* c, doublereal* r) {" << endl; + map out; + m_reactants->writeMultiply("c",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = b->second; + f << " r[" << b->first << "] *= " << rhs << ";" << endl; + } + f << " }" << endl << endl << endl; + } - void ReactionStoichMgr:: - writeMultiplyRevProducts(ostream& f) { - f << " void multiplyRevProducts(const doublereal* c, doublereal* r) {" << endl; - map out; - m_revproducts->writeMultiply("c",out); - map::iterator b; - for (b = out.begin(); b != out.end(); ++b) { - string rhs = b->second; - f << " r[" << b->first << "] *= " << rhs << ";" << endl; - } - f << " }" << endl << endl << endl; - } + void ReactionStoichMgr:: + writeMultiplyRevProducts(ostream& f) { + f << " void multiplyRevProducts(const doublereal* c, doublereal* r) {" << endl; + map out; + m_revproducts->writeMultiply("c",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = b->second; + f << " r[" << b->first << "] *= " << rhs << ";" << endl; + } + f << " }" << endl << endl << endl; + } } diff --git a/Cantera/src/kinetics/ReactionStoichMgr.h b/Cantera/src/kinetics/ReactionStoichMgr.h index d8121d700..203bfc5e4 100644 --- a/Cantera/src/kinetics/ReactionStoichMgr.h +++ b/Cantera/src/kinetics/ReactionStoichMgr.h @@ -5,9 +5,9 @@ */ /* - * $Author: hkmoffa $ - * $Revision: 1.2 $ - * $Date: 2007/06/12 14:20:02 $ + * $Author$ + * $Revision$ + * $Date$ */ #ifndef CT_RXN_STOICH diff --git a/Cantera/src/kinetics/RxnRates.h b/Cantera/src/kinetics/RxnRates.h index a17f549d7..44f0c6a70 100644 --- a/Cantera/src/kinetics/RxnRates.h +++ b/Cantera/src/kinetics/RxnRates.h @@ -3,9 +3,9 @@ * */ -/* $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ +/* $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/StoichManager.h b/Cantera/src/kinetics/StoichManager.h index d699c68d9..56327d85a 100644 --- a/Cantera/src/kinetics/StoichManager.h +++ b/Cantera/src/kinetics/StoichManager.h @@ -1,9 +1,9 @@ /** * @file StoichManager.h * - * $Author: hkmoffa $ - * $Revision: 1.2 $ - * $Date: 2007/05/17 21:17:16 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology @@ -16,808 +16,833 @@ namespace Cantera { - /** - * @defgroup Stoichiometry Stoichiometry - * - * Note: these classes are designed for internal use in class - * ReactionStoichManager. - * - * The classes defined here implement simple operations that are - * used by class ReactionStoichManager to compute things like - * rates of progress, species production rates, etc. In general, a - * reaction mechanism may involve many species and many reactions, - * but any given reaction typically only involves a few species as - * reactants, and a few as products. Therefore, the matrix of - * stoichiometric coefficients is very sparse. Not only is it - * sparse, but the non-zero matrix elements often have the value - * 1, and in many cases no more than three coefficients are - * non-zero for the reactants and/or the products. - * - * For the present purposes, we will consider each direction of a - * reversible reaction to be a separate reaction. We often need to - * compute quantities that can formally be written as a matrix - * product of a stoichiometric coefficient matrix and a vector of - * reaction rates. For example, the species creation rates are - * given by - * \f[ - * \dot C_k = \sum_k \nu^{(p)}_{k,i} R_i - * \f] - * where \f$ \nu^{(p)_{k,i}} \f$ is the product-side stoichiometric - * coefficient of species \a k in reaction \a i. - * This could be done be straightforward matrix multiplication, - * but would be inefficient, since most of the matrix elements - * of \f$ \nu^{(p)}_{k,i} \f$ are zero. We could do better by - * using sparse-matrix algorithms to compute this product. - * - * If the reactions are general ones, with non-integral stoichiometric - * coefficients, this is about as good as we can do. But we are - * particularly concerned here with the performance for very large - * reaction mechanisms, which are usually composed of elementary - * reactions, which have integral stoichiometric - * coefficients. Furthermore, very few elementary reactions involve more - * than 3 product or reactant molecules. This means that instead of + /** + * @defgroup Stoichiometry Stoichiometry + * + * Note: these classes are designed for internal use in class + * ReactionStoichManager. + * + * The classes defined here implement simple operations that are + * used by class ReactionStoichManager to compute things like + * rates of progress, species production rates, etc. In general, a + * reaction mechanism may involve many species and many reactions, + * but any given reaction typically only involves a few species as + * reactants, and a few as products. Therefore, the matrix of + * stoichiometric coefficients is very sparse. Not only is it + * sparse, but the non-zero matrix elements often have the value + * 1, and in many cases no more than three coefficients are + * non-zero for the reactants and/or the products. + * + * For the present purposes, we will consider each direction of a + * reversible reaction to be a separate reaction. We often need to + * compute quantities that can formally be written as a matrix + * product of a stoichiometric coefficient matrix and a vector of + * reaction rates. For example, the species creation rates are + * given by + * + * \f[ + * \dot C_k = \sum_k \nu^{(p)}_{k,i} R_i + * \f] + * + * where \f$ \nu^{(p)_{k,i}} \f$ is the product-side stoichiometric + * coefficient of species \a k in reaction \a i. + * This could be done be straightforward matrix multiplication, + * but would be inefficient, since most of the matrix elements + * of \f$ \nu^{(p)}_{k,i} \f$ are zero. We could do better by + * using sparse-matrix algorithms to compute this product. + * + * If the reactions are general ones, with non-integral stoichiometric + * coefficients, this is about as good as we can do. But we are + * particularly concerned here with the performance for very large + * reaction mechanisms, which are usually composed of elementary + * reactions, which have integral stoichiometric + * coefficients. Furthermore, very few elementary reactions involve more + * than 3 product or reactant molecules. This means that instead of -But we can do even better if we take account of the special structure -of this matrix for elementary reactions. + But we can do even better if we take account of the special structure + of this matrix for elementary reactions. -involve three or fewer product molecules (or reactant molecules). + involve three or fewer product molecules (or reactant molecules). - * To take advantage of this structure, reactions are divided into -These classes are - * designed to take advantage of this sparse structure when - * computing quantities that can be written as matrix multiplies + * To take advantage of this structure, reactions are divided into + * four groups. + These classes are + * designed to take advantage of this sparse structure when + * computing quantities that can be written as matrix multiplies - * They are designed to explicitly unroll loops over species or reactions for - * Operations on reactions that require knowing the reaction - * stoichiometry. - * This module consists of class StoichManager, and - * classes C1, C2, and C3. Classes C1, C2, and C3 handle operations - * involving one, two, or three species, respectively, in a - * reaction. Instances are instantiated with a reaction number, and n - * species numbers (n = 1 for C1, etc.). All three classes have the - * same interface. - * - * These classes are designed for use by StoichManager, and the - * operations implemented are those needed to efficiently compute - * quantities such as rates of progress, species production rates, - * reaction thermochemistry, etc. The compiler will inline these - * methods into the body of the corresponding StoichManager method, - * and so there is no performance penalty (unless inlining is turned - * off). - * - * To describe the methods, consider class C3 and suppose an instance - * is created with reaction number irxn and species numbers k0, k1, - * and k2. - * - * - multiply(in, out) : out[irxn] is multiplied by - * in[k0] * in[k1] * in[k2] - * - * - power(in, out) : out[irxn] is multiplied by - * (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2) - * - * - incrementReaction(in, out) : out[irxn] is incremented by - * in[k0] + in[k1] + in[k2] - * - * - decrementReaction(in, out) : out[irxn] is decremented by - * in[k0] + in[k1] + in[k2] - * - * - incrementSpecies(in, out) : out[k0], out[k1], and out[k2] - * are all incremented by in[irxn] - * - * - decrementSpecies(in, out) : out[k0], out[k1], and out[k2] - * are all decremented by in[irxn] - * - * The function multiply() is usually used when evaluating the - * forward and reverse rates of progress of reactions. - * The rate constants are usually loaded into out[]. Then - * multply() is called to add in the dependence of the - * species concentrations to yield a forward and reverse rop. - * - * The function incrementSpecies() and its cousin decrementSpecies() - * is used to translate from rates of progress to species production - * rates. The vector in[] is preloaed with the rates of progess of - * all reactions. Then incrementSpecies() is called to - * increment the species production vector, out[], with the rates - * of progress. - * - * The functions incrementReaction() and decrementReaction() are - * used to find the standard state equilibrium constant for - * a reaction. Here, output[] is a vector of length - * number of reactions, usually the standard gibbs free energies - * of reaction, while input, usually the standard state - * gibbs free energies of species, is a vector of length number of - * species. - * - * Note the stoichiometric coefficient for a species in a reaction - * is handled by always assuming it is equal to one and then - * treating reactants and products for a reaction separately. - * Bimolecular reactions involving the identical species are - * treated as involving separate species. - * - * @internal This class should be upgraded to include cases where - * real stoichiometric coefficients are used. Shouldn't be that - * hard to do, and they occur in engineering simulations with some - * regularity. - * - */ + * They are designed to explicitly unroll loops over species or reactions for + * Operations on reactions that require knowing the reaction + * stoichiometry. + * This module consists of class StoichManager, and + * classes C1, C2, and C3. Classes C1, C2, and C3 handle operations + * involving one, two, or three species, respectively, in a + * reaction. Instances are instantiated with a reaction number, and n + * species numbers (n = 1 for C1, etc.). All three classes have the + * same interface. + * + * These classes are designed for use by StoichManager, and the + * operations implemented are those needed to efficiently compute + * quantities such as rates of progress, species production rates, + * reaction thermochemistry, etc. The compiler will inline these + * methods into the body of the corresponding StoichManager method, + * and so there is no performance penalty (unless inlining is turned + * off). + * + * To describe the methods, consider class C3 and suppose an instance + * is created with reaction number irxn and species numbers k0, k1, + * and k2. + * + * - multiply(in, out) : out[irxn] is multiplied by + * in[k0] * in[k1] * in[k2] + * + * - power(in, out) : out[irxn] is multiplied by + * (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2) + * + * - incrementReaction(in, out) : out[irxn] is incremented by + * in[k0] + in[k1] + in[k2] + * + * - decrementReaction(in, out) : out[irxn] is decremented by + * in[k0] + in[k1] + in[k2] + * + * - incrementSpecies(in, out) : out[k0], out[k1], and out[k2] + * are all incremented by in[irxn] + * + * - decrementSpecies(in, out) : out[k0], out[k1], and out[k2] + * are all decremented by in[irxn] + * + * The function multiply() is usually used when evaluating the + * forward and reverse rates of progress of reactions. + * The rate constants are usually loaded into out[]. Then + * multply() is called to add in the dependence of the + * species concentrations to yield a forward and reverse rop. + * + * The function incrementSpecies() and its cousin decrementSpecies() + * is used to translate from rates of progress to species production + * rates. The vector in[] is preloaed with the rates of progess of + * all reactions. Then incrementSpecies() is called to + * increment the species production vector, out[], with the rates + * of progress. + * + * The functions incrementReaction() and decrementReaction() are + * used to find the standard state equilibrium constant for + * a reaction. Here, output[] is a vector of length + * number of reactions, usually the standard gibbs free energies + * of reaction, while input, usually the standard state + * gibbs free energies of species, is a vector of length number of + * species. + * + * Note the stoichiometric coefficient for a species in a reaction + * is handled by always assuming it is equal to one and then + * treating reactants and products for a reaction separately. + * Bimolecular reactions involving the identical species are + * treated as involving separate species. + * + * @internal This class should be upgraded to include cases where + * real stoichiometric coefficients are used. Shouldn't be that + * hard to do, and they occur in engineering simulations with some + * regularity. + * + */ - static doublereal ppow(doublereal x, doublereal order) { - if (x > 0.0) - return std::pow(x, order); - else - return 0.0; + static doublereal ppow(doublereal x, doublereal order) { + if (x > 0.0) + return std::pow(x, order); + else + return 0.0; + } + + inline static std::string fmt(std::string r, int n) { return r + "[" + int2str(n) + "]"; } + + + /** + * Handles one species in a reaction. + * @ingroup Stoichiometry + * @internal + */ + class C1 { + + public: + + C1( int rxn = 0, int ic0 = 0) + : m_rxn (rxn), m_ic0 (ic0) {} + + int data(std::vector& ic) { + ic.resize(1); + ic[0] = m_ic0; + return m_rxn; } - inline static std::string fmt(std::string r, int n) { return r + "[" + int2str(n) + "]"; } + void incrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] += R[m_rxn]; + } + void decrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] -= R[m_rxn]; + } + + void multiply(const doublereal* S, doublereal* R) const { + R[m_rxn] *= S[m_ic0]; + } + + void incrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] += S[m_ic0]; + } + + void decrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] -= S[m_ic0]; + } + + int rxnNumber() const { return m_rxn; } + int speciesIndex(int n) const { return m_ic0; } + int nSpecies() { return 1;} + + void writeMultiply(std::string r, std::map& out) { + out[m_rxn] = fmt(r, m_ic0); + } + + void writeIncrementReaction(std::string r, std::map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0); + } + void writeDecrementReaction(std::string r, std::map& out) { + out[m_rxn] += " - "+fmt(r, m_ic0); + } + + void writeIncrementSpecies(std::string r, std::map& out) { + out[m_ic0] += " + "+fmt(r, m_rxn); + } + void writeDecrementSpecies(std::string r, std::map& out) { + out[m_ic0] += " - "+fmt(r, m_rxn); + } + + private: + int m_rxn; + int m_ic0; + }; + + + + /** + * Handles two species in a single reaction. + * @ingroup Stoichiometry + */ + class C2 { + public: + C2( int rxn = 0, int ic0 = 0, int ic1 = 0) + : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1) {} + + int data(std::vector& ic) { + ic.resize(2); + ic[0] = m_ic0; + ic[1] = m_ic1; + return m_rxn; + } + + void incrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] += R[m_rxn]; + S[m_ic1] += R[m_rxn]; + } + + void decrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] -= R[m_rxn]; + S[m_ic1] -= R[m_rxn]; + } + + void multiply(const doublereal* S, doublereal* R) const { + R[m_rxn] *= S[m_ic0] * S[m_ic1]; + } + + void incrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] += S[m_ic0] + S[m_ic1]; + } + + void decrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] -= (S[m_ic0] + S[m_ic1]); + } + + int rxnNumber() const { return m_rxn; } + int speciesIndex(int n) const { return (n == 0 ? m_ic0 : m_ic1); } + int nSpecies() { return 2;} + + void writeMultiply(std::string r, std::map& out) { + out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1); + } + void writeIncrementReaction(std::string r, std::map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1); + } + void writeDecrementReaction(std::string r, std::map& out) { + out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1); + } + + void writeIncrementSpecies(std::string r, std::map& out) { + std::string s = " + "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + } + void writeDecrementSpecies(std::string r, std::map& out) { + std::string s = " - "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + } + + private: /** - * Handles one species in a reaction. - * @ingroup Stoichiometry - * @internal + * Reaction index -> index into the ROP vector */ - class C1 { - - public: - - C1( int rxn = 0, int ic0 = 0) - : m_rxn (rxn), m_ic0 (ic0) {} - - int data(std::vector& ic) { - ic.resize(1); - ic[0] = m_ic0; - return m_rxn; - } - - void incrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] += R[m_rxn]; - } - - void decrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] -= R[m_rxn]; - } - - void multiply(const doublereal* S, doublereal* R) const { - R[m_rxn] *= S[m_ic0]; - } - - void incrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] += S[m_ic0]; - } - - void decrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] -= S[m_ic0]; - } - - int rxnNumber() const { return m_rxn; } - int speciesIndex(int n) const { return m_ic0; } - int nSpecies() { return 1;} - - void writeMultiply(std::string r, std::map& out) { - out[m_rxn] = fmt(r, m_ic0); - } - - void writeIncrementReaction(std::string r, std::map& out) { - out[m_rxn] += " + "+fmt(r, m_ic0); - } - void writeDecrementReaction(std::string r, std::map& out) { - out[m_rxn] += " - "+fmt(r, m_ic0); - } - - void writeIncrementSpecies(std::string r, std::map& out) { - out[m_ic0] += " + "+fmt(r, m_rxn); - } - void writeDecrementSpecies(std::string r, std::map& out) { - out[m_ic0] += " - "+fmt(r, m_rxn); - } - - private: - int m_rxn; - int m_ic0; - }; - - + int m_rxn; /** - * Handles two species in a single reaction. - * @ingroup Stoichiometry + * Species indecise -> index into the species vector for the + * two species. */ - class C2 { - public: - C2( int rxn = 0, int ic0 = 0, int ic1 = 0) - : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1) {} - - int data(std::vector& ic) { - ic.resize(2); - ic[0] = m_ic0; - ic[1] = m_ic1; - return m_rxn; - } - - void incrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] += R[m_rxn]; - S[m_ic1] += R[m_rxn]; - } - - void decrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] -= R[m_rxn]; - S[m_ic1] -= R[m_rxn]; - } - - void multiply(const doublereal* S, doublereal* R) const { - R[m_rxn] *= S[m_ic0] * S[m_ic1]; - } - - void incrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] += S[m_ic0] + S[m_ic1]; - } - - void decrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] -= (S[m_ic0] + S[m_ic1]); - } - - int rxnNumber() const { return m_rxn; } - int speciesIndex(int n) const { return (n == 0 ? m_ic0 : m_ic1); } - int nSpecies() { return 2;} - - void writeMultiply(std::string r, std::map& out) { - out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1); - } - void writeIncrementReaction(std::string r, std::map& out) { - out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1); - } - void writeDecrementReaction(std::string r, std::map& out) { - out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1); - } - - void writeIncrementSpecies(std::string r, std::map& out) { - std::string s = " + "+fmt(r, m_rxn); - out[m_ic0] += s; - out[m_ic1] += s; - } - void writeDecrementSpecies(std::string r, std::map& out) { - std::string s = " - "+fmt(r, m_rxn); - out[m_ic0] += s; - out[m_ic1] += s; - } - - private: - - /** - * Reaction index -> index into the ROP vector - */ - int m_rxn; - - /** - * Species indecise -> index into the species vector for the - * two species. - */ - int m_ic0, m_ic1; - }; + int m_ic0, m_ic1; + }; - /** - * Handles three species in a reaction. - * @ingroup Stoichiometry - */ - class C3 { - public: - C3( int rxn = 0, int ic0 = 0, int ic1 = 0, int ic2 = 0) - : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), m_ic2 (ic2) {} + /** + * Handles three species in a reaction. + * @ingroup Stoichiometry + */ + class C3 { + public: + C3( int rxn = 0, int ic0 = 0, int ic1 = 0, int ic2 = 0) + : m_rxn (rxn), m_ic0 (ic0), m_ic1 (ic1), m_ic2 (ic2) {} - int data(std::vector& ic) { - ic.resize(3); - ic[0] = m_ic0; - ic[1] = m_ic1; - ic[2] = m_ic2; - return m_rxn; - } + int data(std::vector& ic) { + ic.resize(3); + ic[0] = m_ic0; + ic[1] = m_ic1; + ic[2] = m_ic2; + return m_rxn; + } - void incrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] += R[m_rxn]; - S[m_ic1] += R[m_rxn]; - S[m_ic2] += R[m_rxn]; - } + void incrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] += R[m_rxn]; + S[m_ic1] += R[m_rxn]; + S[m_ic2] += R[m_rxn]; + } - void decrementSpecies(const doublereal* R, doublereal* S) const { - S[m_ic0] -= R[m_rxn]; - S[m_ic1] -= R[m_rxn]; - S[m_ic2] -= R[m_rxn]; - } + void decrementSpecies(const doublereal* R, doublereal* S) const { + S[m_ic0] -= R[m_rxn]; + S[m_ic1] -= R[m_rxn]; + S[m_ic2] -= R[m_rxn]; + } - void multiply(const doublereal* S, doublereal* R) const { - R[m_rxn] *= S[m_ic0] * S[m_ic1] * S[m_ic2]; - } + void multiply(const doublereal* S, doublereal* R) const { + R[m_rxn] *= S[m_ic0] * S[m_ic1] * S[m_ic2]; + } - void incrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2]; - } + void incrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2]; + } - void decrementReaction(const doublereal* S, doublereal* R) const { - R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]); - } + void decrementReaction(const doublereal* S, doublereal* R) const { + R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]); + } - int rxnNumber() const { return m_rxn; } - int speciesIndex(int n) const { return (n == 0 ? m_ic0 : (n == 1 ? m_ic1 : m_ic2)); } - int nSpecies() { return 3;} + int rxnNumber() const { return m_rxn; } + int speciesIndex(int n) const { return (n == 0 ? m_ic0 : (n == 1 ? m_ic1 : m_ic2)); } + int nSpecies() { return 3;} - void writeMultiply(std::string r, std::map& out) { - out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1) + " * " + fmt(r, m_ic2); - } - void writeIncrementReaction(std::string r, std::map& out) { - out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1)+" + "+fmt(r, m_ic2); - } - void writeDecrementReaction(std::string r, std::map& out) { - out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1)+" - "+fmt(r, m_ic2); - } - void writeIncrementSpecies(std::string r, std::map& out) { - std::string s = " + "+fmt(r, m_rxn); - out[m_ic0] += s; - out[m_ic1] += s; - out[m_ic2] += s; - } - void writeDecrementSpecies(std::string r, std::map& out) { - std::string s = " - "+fmt(r, m_rxn); - out[m_ic0] += s; - out[m_ic1] += s; - out[m_ic2] += s; - } - private: - int m_rxn, m_ic0, m_ic1, m_ic2; - }; + void writeMultiply(std::string r, std::map& out) { + out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1) + " * " + fmt(r, m_ic2); + } + void writeIncrementReaction(std::string r, std::map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1)+" + "+fmt(r, m_ic2); + } + void writeDecrementReaction(std::string r, std::map& out) { + out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1)+" - "+fmt(r, m_ic2); + } + void writeIncrementSpecies(std::string r, std::map& out) { + std::string s = " + "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + out[m_ic2] += s; + } + void writeDecrementSpecies(std::string r, std::map& out) { + std::string s = " - "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + out[m_ic2] += s; + } + private: + int m_rxn, m_ic0, m_ic1, m_ic2; + }; - /** - * Handles any number of species in a reaction, including fractional - * stoichiometric coefficients, and arbitrary reaction orders. - * @ingroup Stoichiometry - */ - class C_AnyN { - public: - C_AnyN() : m_rxn (-1) {} + /** + * Handles any number of species in a reaction, including fractional + * stoichiometric coefficients, and arbitrary reaction orders. + * @ingroup Stoichiometry + */ + class C_AnyN { + public: + C_AnyN() : m_rxn (-1) {} - C_AnyN( int rxn, const vector_int& ic, const vector_fp& order, + C_AnyN( int rxn, const vector_int& ic, const vector_fp& order, const vector_fp& stoich) - : m_rxn (rxn) { - m_n = ic.size(); - m_ic.resize(m_n); - m_order.resize(m_n); - m_stoich.resize(m_n); - for (int n = 0; n < m_n; n++) { - m_ic[n] = ic[n]; - m_order[n] = order[n]; - m_stoich[n] = stoich[n]; - } - } - - int data(std::vector& ic) { - ic.resize(m_n); - int n; - for (n = 0; n < m_n; n++) ic[n] = m_ic[n]; - return m_rxn; - } - - doublereal order(int n) const {return m_order[n];} - doublereal stoich(int n) const {return m_stoich[n];} - int speciesIndex(int n) const {return m_ic[n];} - - void multiply(const doublereal* input, doublereal* output) const { - for (int n = 0; n < m_n; n++) { - output[m_rxn] *= - ppow(input[m_ic[n]],m_order[n]); - } - } - - void incrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - for (int n = 0; n < m_n; n++) output[m_ic[n]] += m_stoich[n]*x; - } - - void decrementSpecies(const doublereal* input, - doublereal* output) const { - doublereal x = input[m_rxn]; - for (int n = 0; n < m_n; n++) output[m_ic[n]] -= m_stoich[n]*x; - } - - void incrementReaction(const doublereal* input, - doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] - += m_stoich[n]*input[m_ic[n]]; - } - - void decrementReaction(const doublereal* input, - doublereal* output) const { - for (int n = 0; n < m_n; n++) output[m_rxn] - -= m_stoich[n]*input[m_ic[n]]; - } - - void writeMultiply(std::string r, std::map& out) { - int n; - out[m_rxn] = ""; - for (n = 0; n < m_n; n++) { - if (m_order[n] == 1.0) - out[m_rxn] += fmt(r, m_ic[n]); - else - out[m_rxn] += "pow("+fmt(r, m_ic[n])+","+fp2str(m_order[n])+")"; - if (n < m_n-1) - out[m_rxn] += " * "; - } - } - void writeIncrementReaction(std::string r, std::map& out) { - int n; - for (n = 0; n < m_n; n++) { - out[m_rxn] += " + "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]); - } - } - void writeDecrementReaction(std::string r, std::map& out) { - int n; - for (n = 0; n < m_n; n++) { - out[m_rxn] += " - "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]); - } - } - void writeIncrementSpecies(std::string r, std::map& out) { - std::string s = fmt(r, m_rxn); - int n; - for (n = 0; n < m_n; n++) { - out[m_ic[n]] += " + "+fp2str(m_stoich[n]) + "*" + s; - } - } - - void writeDecrementSpecies(std::string r, std::map& out) { - std::string s = fmt(r, m_rxn); - int n; - for (n = 0; n < m_n; n++) { - out[m_ic[n]] += " - "+fp2str(m_stoich[n]) + "*" + s; - } - } - - private: - int m_n, m_rxn; - vector_int m_ic; - vector_fp m_order; - vector_fp m_stoich; - }; - - - template - inline static void _multiply(InputIter begin, InputIter end, - const Vec1& input, Vec2& output) { - for (; begin != end; ++begin) - begin->multiply(input, output); + : m_rxn (rxn) { + m_n = ic.size(); + m_ic.resize(m_n); + m_order.resize(m_n); + m_stoich.resize(m_n); + for (int n = 0; n < m_n; n++) { + m_ic[n] = ic[n]; + m_order[n] = order[n]; + m_stoich[n] = stoich[n]; + } } - template - inline static void _incrementSpecies(InputIter begin, - InputIter end, const Vec1& input, Vec2& output) { - for (; begin != end; ++begin) - begin->incrementSpecies(input, output); + int data(std::vector& ic) { + ic.resize(m_n); + int n; + for (n = 0; n < m_n; n++) ic[n] = m_ic[n]; + return m_rxn; } - template - inline static void _decrementSpecies(InputIter begin, - InputIter end, const Vec1& input, Vec2& output) { - for (; begin != end; ++begin) - begin->decrementSpecies(input, output); + doublereal order(int n) const {return m_order[n];} + doublereal stoich(int n) const {return m_stoich[n];} + int speciesIndex(int n) const {return m_ic[n];} + + void multiply(const doublereal* input, doublereal* output) const { + doublereal oo; + for (int n = 0; n < m_n; n++) { + oo = m_order[n]; + if (oo != 0.0) { + output[m_rxn] *= ppow(input[m_ic[n]], oo); + } + } } - template - inline static void _incrementReactions(InputIter begin, - InputIter end, const Vec1& input, Vec2& output) { - for (; begin != end; ++begin) - begin->incrementReaction(input, output); + void incrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + for (int n = 0; n < m_n; n++) output[m_ic[n]] += m_stoich[n]*x; } - template - inline static void _decrementReactions(InputIter begin, - InputIter end, const Vec1& input, Vec2& output) { - for (; begin != end; ++begin) - begin->decrementReaction(input, output); + void decrementSpecies(const doublereal* input, + doublereal* output) const { + doublereal x = input[m_rxn]; + for (int n = 0; n < m_n; n++) output[m_ic[n]] -= m_stoich[n]*x; } - - template - inline static void _writeIncrementSpecies(InputIter begin, InputIter end, std::string r, - std::map& out) { - for (; begin != end; ++begin) begin->writeIncrementSpecies(r, out); + void incrementReaction(const doublereal* input, + doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] + += m_stoich[n]*input[m_ic[n]]; } - template - inline static void _writeDecrementSpecies(InputIter begin, InputIter end, std::string r, - std::map& out) { - for (; begin != end; ++begin) begin->writeDecrementSpecies(r, out); + void decrementReaction(const doublereal* input, + doublereal* output) const { + for (int n = 0; n < m_n; n++) output[m_rxn] + -= m_stoich[n]*input[m_ic[n]]; } - template - inline static void _writeIncrementReaction(InputIter begin, InputIter end, std::string r, - std::map& out) { - for (; begin != end; ++begin) begin->writeIncrementReaction(r, out); + void writeMultiply(std::string r, std::map& out) { + int n; + out[m_rxn] = ""; + for (n = 0; n < m_n; n++) { + if (m_order[n] == 1.0) + out[m_rxn] += fmt(r, m_ic[n]); + else + out[m_rxn] += "pow("+fmt(r, m_ic[n])+","+fp2str(m_order[n])+")"; + if (n < m_n-1) + out[m_rxn] += " * "; + } + } + void writeIncrementReaction(std::string r, std::map& out) { + int n; + for (n = 0; n < m_n; n++) { + out[m_rxn] += " + "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]); + } + } + void writeDecrementReaction(std::string r, std::map& out) { + int n; + for (n = 0; n < m_n; n++) { + out[m_rxn] += " - "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]); + } + } + void writeIncrementSpecies(std::string r, std::map& out) { + std::string s = fmt(r, m_rxn); + int n; + for (n = 0; n < m_n; n++) { + out[m_ic[n]] += " + "+fp2str(m_stoich[n]) + "*" + s; + } } - template - inline static void _writeDecrementReaction(InputIter begin, InputIter end, std::string r, - std::map& out) { - for (; begin != end; ++begin) begin->writeDecrementReaction(r, out); + void writeDecrementSpecies(std::string r, std::map& out) { + std::string s = fmt(r, m_rxn); + int n; + for (n = 0; n < m_n; n++) { + out[m_ic[n]] += " - "+fp2str(m_stoich[n]) + "*" + s; + } } - template - inline static void _writeMultiply(InputIter begin, InputIter end, std::string r, - std::map& out) { - for (; begin != end; ++begin) begin->writeMultiply(r, out); - } + private: - /* - * This class handles operations involving the stoichiometric - * coefficients on one side of a reaction (reactant or product) for - * a set of reactions comprising a reaction mechanism. This class is - * used by class ReactionStoichMgr, which contains three instances - * of this class (one to handle operations on the reactions, one for - * the products of reversible reactions, and one for the products of - * irreversible reactions). - * - * This class is designed for use with elementary reactions, or at - * least ones with integral stoichiometric coefficients. Let \f$ M(i) \f$ - * be the number of molecules on the product or reactant side of - * reaction number i. - * \f[ - * r_i = \sum_m^{M_i} s_{k_{m,i}} - * \f] - * To understand the operations performed by this class, let - * \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on - * one side (reactant or product) in reaction i. Then \b N is a sparse - * K by I matrix of stoichiometric coefficients. - * - * The following matrix operations may be carried out with a vector - * S of length K, and a vector R of length I: - * - * - \f$ S = S + N R\f$ (incrementSpecies) - * - \f$ S = S - N R\f$ (decrementSpecies) - * - \f$ R = R + N^T S \f$ (incrementReaction) - * - \f$ R = R - N^T S \f$ (deccrementReaction) - * - * The actual implementation, however, does not compute these - * quantities by matrix multiplication. A faster algorithm is used - * that makes use of the fact that the \b integer-valued N matrix is - * very sparse, and the non-zero terms are small positive integers. - * \f[ - * S_k = R_{i1} + \dots + R_{iM} - * \f] - * where M is the number of molecules, and $\f i(m) \f$ is the - * @ingroup Stoichiometry + //! Length of the m_ic vector + /*! + * This is the number of species which have non-zero entries in either the + * reaction order matrix or the stoichiometric order matrix for this reaction. */ - class StoichManagerN { - public: + int m_n; - /** - * Constructor for the StoichManagerN class. - * - * @internal Consider adding defaulted entries here that supply - * the total number of reactions in the mechanism and the total - * number of species in the species list. Then, we could use those - * numbers to provide error checks during the construction of the - * object. Those numbers would also provide some clarity to the - * purpose and utility of this class. - * - * DGG - the problem is that the number of reactions and species - * are not known initially. - */ - StoichManagerN() {} + //! ID of the reaction corresponding to this stoichiometric manager + /*! + * This is used within the interface to select the + */ + int m_rxn; - /** - * Add a single reaction to the list of reactions that this - * stoichiometric manager object handles. - * - * This function is the same as the add() function below. However, - * the order of each species in the power list expression is - * set to one automatically. - */ - void add(int rxn, const vector_int& k) { - vector_fp order(k.size(), 1.0); - vector_fp stoich(k.size(), 1.0); - add(rxn, k, order, stoich); - } - - void add(int rxn, const vector_int& k, const vector_fp& order) { - vector_fp stoich(k.size(), 1.0); - add(rxn, k, order, stoich); - } - - /** - * Add a single reaction to the list of reactions that this - * stoichiometric manager object handles. - * - * @param rxn Reaction index of the current reaction. This is used - * as an index into vectors which have length n_total_rxn. - * @param k This is a vector of integer values specifying the - * species indecises. The length of this vector species - * the number of different species in the description. - * The value of the entries are the species indices. - * These are used as indexes into vectors which have - * length n_total_species. - * @param order This is a vector of the same length as vector k. - * The order is used for the routine power(), which produces - * a power law expression involving the species vector. - * @param stoich This is used to handle fractional stoichiometric coefficients - * on the product side of irreversible reactions. - */ - void add(int rxn, const vector_int& k, const vector_fp& order, - const vector_fp& stoich) { - m_n[rxn] = static_cast(k.size()); - int ns = stoich.size(); - int n; - bool frac = false; - for (n = 0; n < ns; n++) { - if (stoich[n] != 1.0) frac = true; - } - if (frac) { - m_loc[rxn] = static_cast(m_cn_list.size()); - m_cn_list.push_back(C_AnyN(rxn, k, order, stoich)); - } - else { - switch (k.size()) { - case 1: - m_loc[rxn] = static_cast(m_c1_list.size()); - m_c1_list.push_back(C1(rxn, k[0])); - break; - case 2: - m_loc[rxn] = static_cast(m_c2_list.size()); - m_c2_list.push_back(C2(rxn, k[0], k[1])); - break; - case 3: - m_loc[rxn] = static_cast(m_c3_list.size()); - m_c3_list.push_back(C3(rxn, k[0], k[1], k[2])); - break; - default: - m_loc[rxn] = static_cast(m_cn_list.size()); - m_cn_list.push_back(C_AnyN(rxn, k, order, stoich)); - } - } - } - - void multiply(const doublereal* input, doublereal* output) const { - _multiply(m_c1_list.begin(), m_c1_list.end(), input, output); - _multiply(m_c2_list.begin(), m_c2_list.end(), input, output); - _multiply(m_c3_list.begin(), m_c3_list.end(), input, output); - _multiply(m_cn_list.begin(), m_cn_list.end(), input, output); - } - - void incrementSpecies(const doublereal* input, doublereal* output) const { - _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); - _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); - _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); - _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); - } - - void decrementSpecies(const doublereal* input, doublereal* output) const { - _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); - _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); - _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); - _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); - } - - void incrementReactions(const doublereal* input, doublereal* output) const { - _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); - _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); - _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); - _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); - } - - void decrementReactions(const doublereal* input, doublereal* output) const { - _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); - _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); - _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); - _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); - } - - void writeIncrementSpecies(std::string r, std::map& out) { - _writeIncrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out); - _writeIncrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out); - _writeIncrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out); - _writeIncrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out); - } - - void writeDecrementSpecies(std::string r, std::map& out) { - _writeDecrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out); - _writeDecrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out); - _writeDecrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out); - _writeDecrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out); - } - - void writeIncrementReaction(std::string r, std::map& out) { - _writeIncrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out); - _writeIncrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out); - _writeIncrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out); - _writeIncrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out); - } - - void writeDecrementReaction(std::string r, std::map& out) { - _writeDecrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out); - _writeDecrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out); - _writeDecrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out); - _writeDecrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out); - } - - void writeMultiply(std::string r, std::map& out) { - _writeMultiply(m_c1_list.begin(), m_c1_list.end(), r, out); - _writeMultiply(m_c2_list.begin(), m_c2_list.end(), r, out); - _writeMultiply(m_c3_list.begin(), m_c3_list.end(), r, out); - _writeMultiply(m_cn_list.begin(), m_cn_list.end(), r, out); - } + //! Vector of species which are involved with this stoichiometric manager calculations + /*! + * This is an integer list of species which have non-zero entries in either the + * reaction order matrix or the stoichiometric order matrix for this reaction, m_rxn. + * It's used as the index into the arrays m_order[] and m_stoich[]. + */ + vector_int m_ic; + vector_fp m_order; + vector_fp m_stoich; + }; - private: + template + inline static void _multiply(InputIter begin, InputIter end, + const Vec1& input, Vec2& output) { + for (; begin != end; ++begin) + begin->multiply(input, output); + } - std::vector m_c1_list; - std::vector m_c2_list; - std::vector m_c3_list; - std::vector m_cn_list; - /** - * Std::Mapping with the Reaction Number as key and the Number of species - * as the value. - */ - std::map m_n; - /** - * Std::Mapping with the Reaction Number as key and the placement in the - * vector of reactions list( i.e., m_c1_list[]) as key - */ - std::map m_loc; - }; + template + inline static void _incrementSpecies(InputIter begin, + InputIter end, const Vec1& input, Vec2& output) { + for (; begin != end; ++begin) + begin->incrementSpecies(input, output); + } + + template + inline static void _decrementSpecies(InputIter begin, + InputIter end, const Vec1& input, Vec2& output) { + for (; begin != end; ++begin) + begin->decrementSpecies(input, output); + } + + template + inline static void _incrementReactions(InputIter begin, + InputIter end, const Vec1& input, Vec2& output) { + for (; begin != end; ++begin) + begin->incrementReaction(input, output); + } + + template + inline static void _decrementReactions(InputIter begin, + InputIter end, const Vec1& input, Vec2& output) { + for (; begin != end; ++begin) + begin->decrementReaction(input, output); + } + + + template + inline static void _writeIncrementSpecies(InputIter begin, InputIter end, std::string r, + std::map& out) { + for (; begin != end; ++begin) begin->writeIncrementSpecies(r, out); + } + + template + inline static void _writeDecrementSpecies(InputIter begin, InputIter end, std::string r, + std::map& out) { + for (; begin != end; ++begin) begin->writeDecrementSpecies(r, out); + } + + template + inline static void _writeIncrementReaction(InputIter begin, InputIter end, std::string r, + std::map& out) { + for (; begin != end; ++begin) begin->writeIncrementReaction(r, out); + } + + template + inline static void _writeDecrementReaction(InputIter begin, InputIter end, std::string r, + std::map& out) { + for (; begin != end; ++begin) begin->writeDecrementReaction(r, out); + } + + template + inline static void _writeMultiply(InputIter begin, InputIter end, std::string r, + std::map& out) { + for (; begin != end; ++begin) begin->writeMultiply(r, out); + } + + /* + * This class handles operations involving the stoichiometric + * coefficients on one side of a reaction (reactant or product) for + * a set of reactions comprising a reaction mechanism. This class is + * used by class ReactionStoichMgr, which contains three instances + * of this class (one to handle operations on the reactions, one for + * the products of reversible reactions, and one for the products of + * irreversible reactions). + * + * This class is designed for use with elementary reactions, or at + * least ones with integral stoichiometric coefficients. Let \f$ M(i) \f$ + * be the number of molecules on the product or reactant side of + * reaction number i. + * \f[ + * r_i = \sum_m^{M_i} s_{k_{m,i}} + * \f] + * To understand the operations performed by this class, let + * \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on + * one side (reactant or product) in reaction i. Then \b N is a sparse + * K by I matrix of stoichiometric coefficients. + * + * The following matrix operations may be carried out with a vector + * S of length K, and a vector R of length I: + * + * - \f$ S = S + N R\f$ (incrementSpecies) + * - \f$ S = S - N R\f$ (decrementSpecies) + * - \f$ R = R + N^T S \f$ (incrementReaction) + * - \f$ R = R - N^T S \f$ (deccrementReaction) + * + * The actual implementation, however, does not compute these + * quantities by matrix multiplication. A faster algorithm is used + * that makes use of the fact that the \b integer-valued N matrix is + * very sparse, and the non-zero terms are small positive integers. + * \f[ + * S_k = R_{i1} + \dots + R_{iM} + * \f] + * where M is the number of molecules, and $\f i(m) \f$ is the + * @ingroup Stoichiometry + */ + class StoichManagerN { + public: + + /** + * Constructor for the StoichManagerN class. + * + * @internal Consider adding defaulted entries here that supply + * the total number of reactions in the mechanism and the total + * number of species in the species list. Then, we could use those + * numbers to provide error checks during the construction of the + * object. Those numbers would also provide some clarity to the + * purpose and utility of this class. + * + * DGG - the problem is that the number of reactions and species + * are not known initially. + */ + StoichManagerN() {} + + /** + * Add a single reaction to the list of reactions that this + * stoichiometric manager object handles. + * + * This function is the same as the add() function below. However, + * the order of each species in the power list expression is + * set to one automatically. + */ + void add(int rxn, const vector_int& k) { + vector_fp order(k.size(), 1.0); + vector_fp stoich(k.size(), 1.0); + add(rxn, k, order, stoich); + } + + void add(int rxn, const vector_int& k, const vector_fp& order) { + vector_fp stoich(k.size(), 1.0); + add(rxn, k, order, stoich); + } + + + //! Add a single reaction to the list of reactions that this + //! stoichiometric manager object handles. + /*! + * @param rxn Reaction index of the current reaction. This is used + * as an index into vectors which have length n_total_rxn. + * @param k This is a vector of integer values specifying the + * species indecises. The length of this vector species + * the number of different species in the description. + * The value of the entries are the species indices. + * These are used as indexes into vectors which have + * length n_total_species. + * @param order This is a vector of the same length as vector k. + * The order is used for the routine power(), which produces + * a power law expression involving the species vector. + * @param stoich This is used to handle fractional stoichiometric coefficients + * on the product side of irreversible reactions. + */ + void add(int rxn, const vector_int& k, const vector_fp& order, + const vector_fp& stoich) { + m_n[rxn] = static_cast(k.size()); + int ns = stoich.size(); + int n; + bool frac = false; + for (n = 0; n < ns; n++) { + if (stoich[n] != 1.0) frac = true; + } + if (frac) { + m_loc[rxn] = static_cast(m_cn_list.size()); + m_cn_list.push_back(C_AnyN(rxn, k, order, stoich)); + } + else { + switch (k.size()) { + case 1: + m_loc[rxn] = static_cast(m_c1_list.size()); + m_c1_list.push_back(C1(rxn, k[0])); + break; + case 2: + m_loc[rxn] = static_cast(m_c2_list.size()); + m_c2_list.push_back(C2(rxn, k[0], k[1])); + break; + case 3: + m_loc[rxn] = static_cast(m_c3_list.size()); + m_c3_list.push_back(C3(rxn, k[0], k[1], k[2])); + break; + default: + m_loc[rxn] = static_cast(m_cn_list.size()); + m_cn_list.push_back(C_AnyN(rxn, k, order, stoich)); + } + } + } + + void multiply(const doublereal* input, doublereal* output) const { + _multiply(m_c1_list.begin(), m_c1_list.end(), input, output); + _multiply(m_c2_list.begin(), m_c2_list.end(), input, output); + _multiply(m_c3_list.begin(), m_c3_list.end(), input, output); + _multiply(m_cn_list.begin(), m_cn_list.end(), input, output); + } + + void incrementSpecies(const doublereal* input, doublereal* output) const { + _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); + _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); + _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); + _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); + } + + void decrementSpecies(const doublereal* input, doublereal* output) const { + _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); + _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); + _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); + _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); + } + + void incrementReactions(const doublereal* input, doublereal* output) const { + _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); + _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); + _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); + _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); + } + + void decrementReactions(const doublereal* input, doublereal* output) const { + _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); + _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); + _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); + _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); + } + + void writeIncrementSpecies(std::string r, std::map& out) { + _writeIncrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out); + _writeIncrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out); + _writeIncrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out); + _writeIncrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out); + } + + void writeDecrementSpecies(std::string r, std::map& out) { + _writeDecrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out); + _writeDecrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out); + _writeDecrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out); + _writeDecrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out); + } + + void writeIncrementReaction(std::string r, std::map& out) { + _writeIncrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out); + _writeIncrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out); + _writeIncrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out); + _writeIncrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out); + } + + void writeDecrementReaction(std::string r, std::map& out) { + _writeDecrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out); + _writeDecrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out); + _writeDecrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out); + _writeDecrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out); + } + + void writeMultiply(std::string r, std::map& out) { + _writeMultiply(m_c1_list.begin(), m_c1_list.end(), r, out); + _writeMultiply(m_c2_list.begin(), m_c2_list.end(), r, out); + _writeMultiply(m_c3_list.begin(), m_c3_list.end(), r, out); + _writeMultiply(m_cn_list.begin(), m_cn_list.end(), r, out); + } + + + private: + + std::vector m_c1_list; + std::vector m_c2_list; + std::vector m_c3_list; + std::vector m_cn_list; + /** + * Std::Mapping with the Reaction Number as key and the Number of species + * as the value. + */ + std::map m_n; + /** + * Std::Mapping with the Reaction Number as key and the placement in the + * vector of reactions list( i.e., m_c1_list[]) as key + */ + std::map m_loc; + }; #undef INCL_STOICH_WRITER #ifdef INCL_STOICH_WRITER - class StoichWriter { - public: + class StoichWriter { + public: - StoichWriter() {} + StoichWriter() {} - void add(int rxn, const vector_int& k) { - int n, nn = k.size(); - for (n = 0; n < nn; n++) { - if (m_mult[rxn] != "") m_mult[rxn] += " * "; - m_mult[rxn] += "c[" + int2str(k[n]) + "]"; - m_is[k[n]] += " + rop[" + int2str(rxn) + "]"; - m_ds[k[n]] += " - rop[" + int2str(rxn) + "]"; - m_ir[rxn] += " + grt[" + int2str(k[n]) + "]"; - m_dr[rxn] += " - grt[" + int2str(k[n]) + "]"; - } - } + void add(int rxn, const vector_int& k) { + int n, nn = k.size(); + for (n = 0; n < nn; n++) { + if (m_mult[rxn] != "") m_mult[rxn] += " * "; + m_mult[rxn] += "c[" + int2str(k[n]) + "]"; + m_is[k[n]] += " + rop[" + int2str(rxn) + "]"; + m_ds[k[n]] += " - rop[" + int2str(rxn) + "]"; + m_ir[rxn] += " + grt[" + int2str(k[n]) + "]"; + m_dr[rxn] += " - grt[" + int2str(k[n]) + "]"; + } + } - void add(int rxn, const vector_int& k, const vector_fp& order, - const vector_fp& stoich) { - int n, nn = k.size(); - std::string s; - for (n = 0; n < nn; n++) { - if (order[n] == 1.0) - m_mult[rxn] += "*c[" + int2str(k[n]) + "]"; - else - m_mult[rxn] += "*pow(c[" _ int2str(k[n]) + "],"+fp2str(order[n])+")"; - if (stoich[n] == 1.0) { - m_is[k[n]] += " + r[" + int2str(rxn) + "]"; - m_ds[k[n]] += " - r[" + int2str(rxn) + "]"; - m_ir[rxn] += " + g[" + int2str(k[n]) + "]"; - m_dr[rxn] += " - g[" + int2str(k[n]) + "]"; - } - else { - s = fp2str(stoich[n]); - m_is[k[n]] += " + "+s+"*r[" + int2str(rxn) + "]"; - m_ds[k[n]] += " - "+s+"*r[" + int2str(rxn) + "]"; - m_ir[rxn] += " + "+s+"*g[" + int2str(k[n]) + "]"; - m_dr[rxn] += " - "+s+"*g[" + int2str(k[n]) + "]"; - } - } - } + void add(int rxn, const vector_int& k, const vector_fp& order, + const vector_fp& stoich) { + int n, nn = k.size(); + std::string s; + for (n = 0; n < nn; n++) { + if (order[n] == 1.0) + m_mult[rxn] += "*c[" + int2str(k[n]) + "]"; + else + m_mult[rxn] += "*pow(c[" _ int2str(k[n]) + "],"+fp2str(order[n])+")"; + if (stoich[n] == 1.0) { + m_is[k[n]] += " + r[" + int2str(rxn) + "]"; + m_ds[k[n]] += " - r[" + int2str(rxn) + "]"; + m_ir[rxn] += " + g[" + int2str(k[n]) + "]"; + m_dr[rxn] += " - g[" + int2str(k[n]) + "]"; + } + else { + s = fp2str(stoich[n]); + m_is[k[n]] += " + "+s+"*r[" + int2str(rxn) + "]"; + m_ds[k[n]] += " - "+s+"*r[" + int2str(rxn) + "]"; + m_ir[rxn] += " + "+s+"*g[" + int2str(k[n]) + "]"; + m_dr[rxn] += " - "+s+"*g[" + int2str(k[n]) + "]"; + } + } + } - std::string mult(int rxn) { return m_mult[rxn]; } - std::string incrSpec(int k, std::string) { return m_is[k]; } - std::string decrSpec(int k) { return m_ds[k]; } - std::string incrRxn(int rxn) { return m_ir[rxn]; } - std::string decrRxn(int rxn) { return m_dr[rxn]; } + std::string mult(int rxn) { return m_mult[rxn]; } + std::string incrSpec(int k, std::string) { return m_is[k]; } + std::string decrSpec(int k) { return m_ds[k]; } + std::string incrRxn(int rxn) { return m_ir[rxn]; } + std::string decrRxn(int rxn) { return m_dr[rxn]; } - private: - std::map m_mult, m_ir, m_dr, m_is, m_ds; - }; + private: + std::map m_mult, m_ir, m_dr, m_is, m_ds; + }; #endif diff --git a/Cantera/src/kinetics/ThirdBodyMgr.h b/Cantera/src/kinetics/ThirdBodyMgr.h index 63ccea4f8..33ab92c22 100644 --- a/Cantera/src/kinetics/ThirdBodyMgr.h +++ b/Cantera/src/kinetics/ThirdBodyMgr.h @@ -1,9 +1,9 @@ /** * @file ThirdBodyMgr.h * - * $Author: dggoodwin $ - * $Revision: 1.1 $ - * $Date: 2007/05/04 14:27:23 $ + * $Author$ + * $Revision$ + * $Date$ */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/importKinetics.cpp b/Cantera/src/kinetics/importKinetics.cpp index 7641bdd9b..fd6b5a1c4 100644 --- a/Cantera/src/kinetics/importKinetics.cpp +++ b/Cantera/src/kinetics/importKinetics.cpp @@ -10,9 +10,9 @@ * from the ctml tree structures. */ -/* $Author: hkmoffa $ - * $Revision: 1.10 $ - * $Date: 2009/03/13 03:23:20 $ +/* $Author$ + * $Revision$ + * $Date$ */ // Copyright 2002 California Institute of Technology @@ -50,32 +50,32 @@ using namespace std; namespace Cantera { -//! these are all used to check for duplicate reactions -class rxninfo { -public: - //! rdata - std::vector< std::map > m_rdata; - //! string name - std::vector m_eqn; - //! string vector of ints - std::vector m_dup; - //! string vector of ints - std::vector m_nr; - //! string vector of ints - std::vector m_typ; - //! vector of bools. - std::vector m_rev; - ~rxninfo() { - m_eqn.clear(); - m_dup.clear(); - m_nr.clear(); - m_typ.clear(); - m_rdata.clear(); - } - bool installReaction(int i, const XML_Node& r, Kinetics* k, - std::string default_phase, int rule, - bool validate_rxn) ; -}; + //! these are all used to check for duplicate reactions + class rxninfo { + public: + //! rdata + std::vector< std::map > m_rdata; + //! string name + std::vector m_eqn; + //! string vector of ints + std::vector m_dup; + //! string vector of ints + std::vector m_nr; + //! string vector of ints + std::vector m_typ; + //! vector of bools. + std::vector m_rev; + ~rxninfo() { + m_eqn.clear(); + m_dup.clear(); + m_nr.clear(); + m_typ.clear(); + m_rdata.clear(); + } + bool installReaction(int i, const XML_Node& r, Kinetics* k, + std::string default_phase, int rule, + bool validate_rxn) ; + }; @@ -266,29 +266,29 @@ public: * Check to see if reaction orders have been specified. */ if (rp == 1 && rxn.hasChild("order")) { - vector ord; - rxn.getChildren("order",ord); - int norder = static_cast(ord.size()); - int loc; - doublereal forder; - for (int nn = 0; nn < norder; nn++) { - const XML_Node& oo = *ord[nn]; - string sp = oo["species"]; - loc = speciesMap[sp]; - if (loc == 0) - throw CanteraError("getReagents", - "reaction order specified for non-reactantt: " - +sp); - forder = fpValue(oo()); - if (forder < 0.0) { - throw CanteraError("getReagents", - "reaction order must be non-negative"); - } - // replace the stoichiometric coefficient - // stored above in 'order' with the specified - // reaction order - order[loc-1] = forder; - } + vector ord; + rxn.getChildren("order",ord); + int norder = static_cast(ord.size()); + int loc; + doublereal forder; + for (int nn = 0; nn < norder; nn++) { + const XML_Node& oo = *ord[nn]; + string sp = oo["species"]; + loc = speciesMap[sp]; + if (loc == 0) + throw CanteraError("getReagents", + "reaction order specified for non-reactantt: " + +sp); + forder = fpValue(oo()); + if (forder < 0.0) { + throw CanteraError("getReagents", + "reaction order must be non-negative"); + } + // replace the stoichiometric coefficient + // stored above in 'order' with the specified + // reaction order + order[loc-1] = forder; + } } return true; } @@ -314,117 +314,117 @@ public: E /= GasConstant; } - /** - * getStick() processes the XML element called Stick that specifies - * the sticking coefficient reaction. This routine will - * translate the sticking coefficient value into a "normal" - * rate constant for the surface reaction. - * - * Output - * ----------- - * Output is the normal Arrhenius expressions for a surface - * reaction rate constant. - * - * A - units such that rate of rxn has kmol/m^2/s when - * A is multiplied by activity concentrations of - * reactants in the normal manner. - * n - unitless - * E - Units 1/Kelvin + /** + * getStick() processes the XML element called Stick that specifies + * the sticking coefficient reaction. This routine will + * translate the sticking coefficient value into a "normal" + * rate constant for the surface reaction. + * + * Output + * ----------- + * Output is the normal Arrhenius expressions for a surface + * reaction rate constant. + * + * A - units such that rate of rxn has kmol/m^2/s when + * A is multiplied by activity concentrations of + * reactants in the normal manner. + * n - unitless + * E - Units 1/Kelvin + */ + static void getStick(const XML_Node& node, Kinetics& kin, + ReactionData& r, doublereal& A, doublereal& b, doublereal& E) { + int nr = r.reactants.size(); + int k, klocal, not_surf = 0; + int np = 0; + doublereal f = 1.0; + doublereal order; + /* + * species is the name of the special reactant whose surface + * flux rate will be calculated. + * isp = species # in the local phase + * ispKinetics = species # in the kinetics object + * ispPhaseIndex = phase # of the special species */ - static void getStick(const XML_Node& node, Kinetics& kin, - ReactionData& r, doublereal& A, doublereal& b, doublereal& E) { - int nr = r.reactants.size(); - int k, klocal, not_surf = 0; - int np = 0; - doublereal f = 1.0; - doublereal order; - /* - * species is the name of the special reactant whose surface - * flux rate will be calculated. - * isp = species # in the local phase - * ispKinetics = species # in the kinetics object - * ispPhaseIndex = phase # of the special species - */ - string spname = node["species"]; - ThermoPhase& th = kin.speciesPhase(spname); - int isp = th.speciesIndex(spname); - int ispKinetics = kin.kineticsSpeciesIndex(spname); - int ispPhaseIndex = kin.speciesPhaseIndex(ispKinetics); + string spname = node["species"]; + ThermoPhase& th = kin.speciesPhase(spname); + int isp = th.speciesIndex(spname); + int ispKinetics = kin.kineticsSpeciesIndex(spname); + int ispPhaseIndex = kin.speciesPhaseIndex(ispKinetics); - double ispMW = th.molecularWeights()[isp]; - double sc; + double ispMW = th.molecularWeights()[isp]; + double sc; - // loop over the reactants - for (int n = 0; n < nr; n++) { - k = r.reactants[n]; - order = r.order[n]; // stoich coeff + // loop over the reactants + for (int n = 0; n < nr; n++) { + k = r.reactants[n]; + order = r.rorder[n]; // stoich coeff - // get the phase species k belongs to - np = kin.speciesPhaseIndex(k); - const ThermoPhase& p = kin.thermo(np); + // get the phase species k belongs to + np = kin.speciesPhaseIndex(k); + const ThermoPhase& p = kin.thermo(np); - // get the local index of species k in this phase - klocal = p.speciesIndex(kin.kineticsSpeciesName(k)); + // get the local index of species k in this phase + klocal = p.speciesIndex(kin.kineticsSpeciesName(k)); - // if it is a surface species, divide f by the standard - // concentration for this species, in order to convert - // from concentration units used in the law of mass action - // to coverages used in the sticking probability - // expression - if (p.eosType() == cSurf || p.eosType() == cEdge) { - sc = p.standardConcentration(klocal); - f /= pow(sc, order); - } - // Otherwise: - else { - // We only allow one species to be in the phase - // containing the special sticking coefficient - // species. - if (ispPhaseIndex == np) { - not_surf++; - } - // Other bulk phase species on the other side - // of ther interface are treated like surface - // species. - else { - sc = p.standardConcentration(klocal); - f /= pow(sc, order); - } - } - } - if (not_surf != 1) { - throw CanteraError("getStick", - "reaction probabilities can only be used in " - "reactions with exactly 1 gas/liquid species."); - } - - doublereal cbar = sqrt(8.0*GasConstant/(Pi*ispMW)); - A = 0.25 * getFloat(node, "A", "toSI") * cbar * f; - b = getFloat(node, "b") + 0.5; - E = getFloat(node, "E", "actEnergy"); - E /= GasConstant; - } - - static void getCoverageDependence(const node_t& node, - thermo_t& surfphase, ReactionData& rdata) { - vector cov; - node.getChildren("coverage", cov); - int k, nc = static_cast(cov.size()); - doublereal e; - string spname; - if (nc > 0) { - for (int n = 0; n < nc; n++) { - const XML_Node& cnode = *cov[n]; - spname = cnode["species"]; - k = surfphase.speciesIndex(spname); - rdata.cov.push_back(doublereal(k)); - rdata.cov.push_back(getFloat(cnode, "a")); - rdata.cov.push_back(getFloat(cnode, "m")); - e = getFloat(cnode, "e", "actEnergy"); - rdata.cov.push_back(e/GasConstant); - } - } + // if it is a surface species, divide f by the standard + // concentration for this species, in order to convert + // from concentration units used in the law of mass action + // to coverages used in the sticking probability + // expression + if (p.eosType() == cSurf || p.eosType() == cEdge) { + sc = p.standardConcentration(klocal); + f /= pow(sc, order); + } + // Otherwise: + else { + // We only allow one species to be in the phase + // containing the special sticking coefficient + // species. + if (ispPhaseIndex == np) { + not_surf++; + } + // Other bulk phase species on the other side + // of ther interface are treated like surface + // species. + else { + sc = p.standardConcentration(klocal); + f /= pow(sc, order); + } + } } + if (not_surf != 1) { + throw CanteraError("getStick", + "reaction probabilities can only be used in " + "reactions with exactly 1 gas/liquid species."); + } + + doublereal cbar = sqrt(8.0*GasConstant/(Pi*ispMW)); + A = 0.25 * getFloat(node, "A", "toSI") * cbar * f; + b = getFloat(node, "b") + 0.5; + E = getFloat(node, "E", "actEnergy"); + E /= GasConstant; + } + + static void getCoverageDependence(const node_t& node, + thermo_t& surfphase, ReactionData& rdata) { + vector cov; + node.getChildren("coverage", cov); + int k, nc = static_cast(cov.size()); + doublereal e; + string spname; + if (nc > 0) { + for (int n = 0; n < nc; n++) { + const XML_Node& cnode = *cov[n]; + spname = cnode["species"]; + k = surfphase.speciesIndex(spname); + rdata.cov.push_back(doublereal(k)); + rdata.cov.push_back(getFloat(cnode, "a")); + rdata.cov.push_back(getFloat(cnode, "m")); + e = getFloat(cnode, "e", "actEnergy"); + rdata.cov.push_back(e/GasConstant); + } + } + } //! Get falloff parameters for a reaction. @@ -434,9 +434,9 @@ public: * * * @verbatim - 0.5 73.2 5000. 9999. - @endverbatim - */ + 0.5 73.2 5000. 9999. + @endverbatim + */ static void getFalloff(const node_t& f, ReactionData& rdata) { string type = f["type"]; vector p; @@ -469,7 +469,7 @@ public: } else { throw CanteraError("getFalloff()", "Troe parameterization is specified by number of pararameters, " - + int2str(np) + ", is not equal to 3 or 4"); + + int2str(np) + ", is not equal to 3 or 4"); } } else if (type == "SRI") { if (np == 5) { @@ -487,7 +487,7 @@ public: } } else { throw CanteraError("getFalloff()", "SRI parameterization is specified by number of pararameters, " - + int2str(np) + ", is not equal to 3 or 5"); + + int2str(np) + ", is not equal to 3 or 5"); } } rdata.falloffParameters = c; @@ -645,8 +645,8 @@ public: * @ingroup kineticsmgr */ bool rxninfo::installReaction(int i, const XML_Node& r, Kinetics* k, - string default_phase, int rule, - bool validate_rxn) { + string default_phase, int rule, + bool validate_rxn) { Kinetics& kin = *k; @@ -701,7 +701,7 @@ public: // get the reactants bool ok = getReagents(r, kin, 1, default_phase, rdata.reactants, - rdata.rstoich, rdata.order, rule); + rdata.rstoich, rdata.rorder, rule); //cout << "Reactants: " << endl; //int npp = rdata.reactants.size(); //int nj; @@ -712,9 +712,8 @@ public: /* * Get the products. We store the id of products in rdata.products */ - vector_fp dummy; ok = ok && getReagents(r, kin, -1, default_phase, rdata.products, - rdata.pstoich, dummy, rule); + rdata.pstoich, rdata.porder, rule); //cout << "Products: " << endl;npp = rdata.products.size(); //for (nj = 0; nj < npp; nj++) { // cout << rdata.products[nj] << " " << rdata.pstoich[nj] << endl; @@ -749,6 +748,55 @@ public: rdata.global = true; } + /* + * Some reactions can be elementary reactions but have fractional + * stoichiometries wrt to some products and reactants. An + * example of these are solid reactions involving phase transformations. + * Species with fractional stoichiometries must be from single-species + * phases with unity activities. For these reactions set + * the bool isReversibleWithFrac to true. + */ + if (rdata.reversible == true) { + int np = rdata.products.size(); + for (int i = 0; i < np; i++) { + int k = rdata.products[i]; + double po = rdata.porder[i]; + AssertTrace(po == rdata.pstoich[i]); + double chk = po - 1.0 * int(po); + if (chk != 0.0) { + /* + * put in a check here that k is a single species phase. + */ + thermo_t &thref = kin.speciesPhase(k); + if (thref.nSpecies() == 1) { + rdata.porder[i] = 0.0; + } + + rdata.isReversibleWithFrac = true; + + } + } + int nr = rdata.reactants.size(); + for (int i = 0; i < nr; i++) { + int k = rdata.reactants[i]; + double ro = rdata.rorder[i]; + AssertTrace(ro == rdata.rstoich[i]); + double chk = ro - 1.0 * int(ro); + if (chk != 0.0) { + /* + * put in a check here that k is a single species phase. + */ + thermo_t &thref = kin.speciesPhase(k); + if (thref.nSpecies() == 1) { + rdata.rorder[i] = 0.0; + } + + rdata.isReversibleWithFrac = true; + + } + } + } + /* * Search the reaction element for the attribute "type". * If found, then branch on the type, to fill in appropriate @@ -934,7 +982,7 @@ public: const XML_Node* r = allrxns[i]; if (r) { if (_rxns->installReaction(itot, *r, &kin, - default_phase, rxnrule, check_for_duplicates)) ++itot; + default_phase, rxnrule, check_for_duplicates)) ++itot; } } } @@ -968,7 +1016,7 @@ public: */ if ((rxid >= imin) && (rxid <= imax)) { if (_rxns->installReaction(itot, *r, &kin, - default_phase, rxnrule, check_for_duplicates)) ++itot; + default_phase, rxnrule, check_for_duplicates)) ++itot; } } } diff --git a/Cantera/src/kinetics/importKinetics.h b/Cantera/src/kinetics/importKinetics.h index 8735ac18d..404937bc7 100644 --- a/Cantera/src/kinetics/importKinetics.h +++ b/Cantera/src/kinetics/importKinetics.h @@ -10,9 +10,9 @@ * from the ctml tree structures. */ /* - * $Author: hkmoffa $ - * $Revision: 1.2 $ - * $Date: 2007/06/04 23:09:19 $ + * $Author$ + * $Revision$ + * $Date$ * */ diff --git a/Cantera/src/kinetics/reaction_defs.h b/Cantera/src/kinetics/reaction_defs.h index 8e6123915..eb36dc544 100644 --- a/Cantera/src/kinetics/reaction_defs.h +++ b/Cantera/src/kinetics/reaction_defs.h @@ -4,9 +4,9 @@ */ /* - * $Author: dggoodwin $ - * $Date: 2007/05/04 14:27:24 $ - * $Revision: 1.1 $ + * $Author$ + * $Date$ + * $Revision$ * */ diff --git a/Cantera/src/kinetics/solveSP.cpp b/Cantera/src/kinetics/solveSP.cpp index a91f0a3c9..a44757365 100644 --- a/Cantera/src/kinetics/solveSP.cpp +++ b/Cantera/src/kinetics/solveSP.cpp @@ -2,7 +2,7 @@ * @file: solveSP.cpp Implicit surface site concentration solver */ /* - * $Id: solveSP.cpp,v 1.4 2008/12/17 17:09:37 hkmoffa Exp $ + * $Id$ */ /* * Copywrite 2004 Sandia Corporation. Under the terms of Contract diff --git a/Cantera/src/kinetics/solveSP.h b/Cantera/src/kinetics/solveSP.h index 7caf9b334..f4d679d3e 100644 --- a/Cantera/src/kinetics/solveSP.h +++ b/Cantera/src/kinetics/solveSP.h @@ -4,7 +4,7 @@ * (see \ref chemkinetics and class \link Cantera::solveSP solveSP\endlink). */ /* - * $Id: solveSP.h,v 1.3 2007/12/30 04:19:38 dggoodwin Exp $ + * $Id$ */ /* * Copywrite 2004 Sandia Corporation. Under the terms of Contract