From b0715adf56336f417e959bdb0a07cc516bbc51e0 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 18 Aug 2003 17:54:11 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/src/GasKinetics.cpp | 30 +++++----- Cantera/src/InterfaceKinetics.cpp | 57 +++++++++++++++---- Cantera/src/InterfaceKinetics.h | 9 ++- Cantera/src/Makefile.in | 2 +- Cantera/src/ReactionData.h | 1 + Cantera/src/ReactionStoichMgr.h | 86 ++++++++++++++++++++++++++++- Cantera/src/RxnRates.h | 84 +++++++++++++++++++++++++++- Cantera/src/StoichManager.h | 5 +- Cantera/src/converters/CKReader.cpp | 2 +- Cantera/src/converters/ck2ct.cpp | 4 +- Cantera/src/importCTML.cpp | 5 -- Cantera/src/importCTML.h | 14 +++-- Cantera/src/reaction_defs.h | 3 +- 13 files changed, 253 insertions(+), 49 deletions(-) diff --git a/Cantera/src/GasKinetics.cpp b/Cantera/src/GasKinetics.cpp index c258774e9..ba263552d 100755 --- a/Cantera/src/GasKinetics.cpp +++ b/Cantera/src/GasKinetics.cpp @@ -68,7 +68,7 @@ namespace Cantera { _update_rates_T() { doublereal T = thermo().temperature(); m_kdata->m_logc0 = log(thermo().standardConcentration()); - if (fabs(T - m_kdata->m_temp) > m_dt_threshold) { + if (fabs(T - m_kdata->m_temp) > 0.0) { // m_dt_threshold) { doublereal logT = log(T); //m_kdata->m_logp0 - logT; m_rates.update(T, logT, m_kdata->m_rfn.begin()); @@ -79,20 +79,20 @@ namespace Cantera { updateKc(); m_kdata->m_ROP_ok = false; } - else { - doublereal logT = log(T); - doublereal dT = T - m_kdata->m_temp; - //m_kdata->m_logc0 = m_kdata->m_logp0 - logT; - m_rates.update_dT(T, logT, dT, m_kdata->m_rfn.begin()); - m_falloff_low_rates.update_dT(T, logT, dT, - m_kdata->m_rfn_low.begin()); - m_falloff_high_rates.update_dT(T, logT, dT, - m_kdata->m_rfn_high.begin()); - m_falloffn.updateTemp(T, m_kdata->falloff_work.begin()); - m_kdata->m_temp = T; - updateKc(); - m_kdata->m_ROP_ok = false; - } +// else { +// doublereal logT = log(T); +// doublereal dT = T - m_kdata->m_temp; +// //m_kdata->m_logc0 = m_kdata->m_logp0 - logT; +// m_rates.update_dT(T, logT, dT, m_kdata->m_rfn.begin()); +// m_falloff_low_rates.update_dT(T, logT, dT, +// m_kdata->m_rfn_low.begin()); +// m_falloff_high_rates.update_dT(T, logT, dT, +// m_kdata->m_rfn_high.begin()); +// m_falloffn.updateTemp(T, m_kdata->falloff_work.begin()); +// m_kdata->m_temp = T; +// updateKc(); +// m_kdata->m_ROP_ok = false; +// } }; diff --git a/Cantera/src/InterfaceKinetics.cpp b/Cantera/src/InterfaceKinetics.cpp index 3c0f99e0a..6b4b038d9 100644 --- a/Cantera/src/InterfaceKinetics.cpp +++ b/Cantera/src/InterfaceKinetics.cpp @@ -194,6 +194,7 @@ namespace Cantera { m_redo_rates(false), m_nirrev(0), m_nrev(0), + m_surf(0), m_integrator(0), m_finalized(false) { @@ -218,10 +219,10 @@ namespace Cantera { void InterfaceKinetics:: _update_rates_T() { _update_rates_phi(); - doublereal T = thermo().temperature(); + doublereal T = thermo(surfacePhaseIndex()).temperature(); if (T != m_kdata->m_temp || m_redo_rates) { - doublereal logT = log(T); - m_rates.update(T, logT, m_kdata->m_rfn.begin()); + m_kdata->m_logtemp = log(T); + m_rates.update(T, m_kdata->m_logtemp, m_kdata->m_rfn.begin()); applyButlerVolmerCorrection(m_kdata->m_rfn.begin()); m_kdata->m_temp = T; updateKc(); @@ -252,6 +253,16 @@ namespace Cantera { void InterfaceKinetics:: _update_rates_C() { int n; + + /** + * First evaluate the coverage-dependent terms in the reaction + * rates. + */ + m_surf->getCoverages(m_conc.begin()); + m_rates.update_C(m_conc.begin()); + m_rates.update(m_kdata->m_temp, + m_kdata->m_logtemp, m_kdata->m_rfn.begin()); + int np = nPhases(); for (n = 0; n < np; n++) { /* @@ -452,12 +463,24 @@ namespace Cantera { void InterfaceKinetics:: addReaction(const ReactionData& r) { - if (r.reactionType == ELEMENTARY_RXN) - addElementaryReaction(r); - if (r.reactionType == SURFACE_RXN) - addElementaryReaction(r); - else if (r.reactionType == GLOBAL_RXN) + int nr = r.reactants.size(); + bool isglobal = false; + for (int n = 0; n < nr; n++) { + if (r.rstoich[n] != int(r.order[n])) { + isglobal = true; break; + } + } + if (isglobal) addGlobalReaction(r); + else + addElementaryReaction(r); + +// if (r.reactionType == ELEMENTARY_RXN) +// addElementaryReaction(r); +// if (r.reactionType == SURFACE_RXN) +// addElementaryReaction(r); +// else if (r.reactionType == GLOBAL_RXN) +// addGlobalReaction(r); // operations common to all reaction types installReagents( r ); @@ -471,9 +494,12 @@ namespace Cantera { addElementaryReaction(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, r.rateCoeffParameters.size(), - r.rateCoeffParameters.begin() ); + r.rateCoeffType, rp.size(), + rp.begin() ); // store activation energy m_E.push_back(r.rateCoeffParameters[2]); // add constant term to rate coeff value vector @@ -487,9 +513,12 @@ namespace Cantera { 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, r.rateCoeffParameters.size(), - r.rateCoeffParameters.begin() ); + r.rateCoeffType, rp.size(), + rp.begin() ); // add constant term to rate coeff value vector m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); @@ -597,6 +626,10 @@ namespace Cantera { */ void InterfaceKinetics::finalize() { m_rwork.resize(nReactions()); + int ks = surfacePhaseIndex(); + if (ks < 0) throw CanteraError("InterfaceKinetics::finalize", + "no surface phase is present."); + m_surf = (SurfPhase*)&thermo(ks); m_finalized = true; } diff --git a/Cantera/src/InterfaceKinetics.h b/Cantera/src/InterfaceKinetics.h index 392de81ce..4f528db16 100644 --- a/Cantera/src/InterfaceKinetics.h +++ b/Cantera/src/InterfaceKinetics.h @@ -31,6 +31,7 @@ namespace Cantera { class ReactionData; class InterfaceKineticsData; class ThermoPhase; + class SurfPhase; class ImplicitSurfChem; /** @@ -40,7 +41,7 @@ namespace Cantera { public: InterfaceKineticsData() : m_ROP_ok(false), - m_temp(0.0) + m_temp(0.0), m_logtemp(0.0) {} virtual ~InterfaceKineticsData(){} @@ -49,7 +50,7 @@ namespace Cantera { array_fp m_rfn_low, m_rfn_high; bool m_ROP_ok; - doublereal m_temp; + doublereal m_temp, m_logtemp; vector_fp m_rfn; vector_fp m_rkcn; }; @@ -171,7 +172,8 @@ namespace Cantera { int m_kk; - Rate1 m_rates; + Rate1 m_rates; + //Rate1 m_rates; bool m_redo_rates; mutable map > m_index; @@ -221,6 +223,7 @@ namespace Cantera { vector_fp m_rwork; vector_fp m_E; + SurfPhase* m_surf; ImplicitSurfChem* m_integrator; private: diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index fe2dcaca8..9ef79f5fe 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -29,7 +29,7 @@ BASE = State.o Elements.o Constituents.o stringUtils.o misc.o importCTML.o THERMO = $(BASE) ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o SolidCompound.o SpeciesThermoFactory.o ThermoFactory.o # homogeneous kinetics -KINETICS = GRI_30_Kinetics.o KineticsFactory.o GasKinetics.o FalloffFactory.o GasKineticsWriter.o \ +KINETICS = GRI_30_Kinetics.o KineticsFactory.o GasKinetics.o FalloffFactory.o \ ReactionStoichMgr.o $(THERMO) # heterogeneous kinetics diff --git a/Cantera/src/ReactionData.h b/Cantera/src/ReactionData.h index 0b4bb5944..7c2d8d0cb 100755 --- a/Cantera/src/ReactionData.h +++ b/Cantera/src/ReactionData.h @@ -55,6 +55,7 @@ namespace Cantera { int error; string equation; doublereal default_3b_eff; + vector_fp cov; }; } diff --git a/Cantera/src/ReactionStoichMgr.h b/Cantera/src/ReactionStoichMgr.h index 7925d4207..a2bece857 100644 --- a/Cantera/src/ReactionStoichMgr.h +++ b/Cantera/src/ReactionStoichMgr.h @@ -1,12 +1,24 @@ +/** + * @file ReactionStoichMgr.h + * + * $Author$ + * $Revision$ + * $Date$ + */ + #ifndef CT_RXN_STOICH #define CT_RXN_STOICH + #include "ct_defs.h" namespace Cantera { class StoichManagerN; + /** + * This class handles calculations involving reaction stoichiometry. + */ class ReactionStoichMgr { public: @@ -14,18 +26,88 @@ namespace Cantera { ReactionStoichMgr(); virtual ~ReactionStoichMgr(); + /** + * Add a reaction with specified, possibly non-integral, reaction orders. + * @param rxn Reaction number + * @param reactants vector of integer reactant indices + * @param reactants vector of integer product indices + * @param reversible true if the reaction is reversible, false otherwise + * @param fwdOrder reaction orders for the reactants. This vector must + * be the same length as 'reactants,' and the reaction orders are for the + * species with index in the corresponding location in 'reactants.' + */ void add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible, const vector_fp& fwdOrder); + bool reversible, const vector_fp& fwdOrder); + + /** + * Add a reaction with mass-action kinetics. + * @param rxn Reaction number + * @param reactants vector of integer reactant indices + * @param reactants vector of integer product indices + * @param reversible true if the reaction is reversible, false otherwise + */ void add(int rxn, const vector_int& reactants, const vector_int& products, bool reversible); + /** + * Given the arrays of the forward and reverse rates of progress for all reactions, + * compute the species creation rates and return them in array c. + */ void getCreationRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c); + + /** + * Given the arrays of the forward and reverse rates of progress for all reactions, + * compute the species destruction rates and return them in array d. + */ void getDestructionRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* d); + + /** + * Given the array of the net rates of progress for all reactions, + * compute the species net production rates and return them in array w. + */ void getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w); + + /** + * Given an array of species properties 'g', return in array 'dg' the change in this quantity + * in the reactions. Array 'g' must have a length at least as great + * as the number of species, and array 'dg' must have a length + * as great as the total number of reactions. + */ void getReactionDelta(int nr, const doublereal* g, doublereal* dg); + + /** + * Given an array of species properties 'g', return in array + * 'dg' the change in this quantity in the reversible + * reactions. Array 'g' must have a length at least as great + * as the number of species, and array 'dg' must have a length + * as great as the total number of reactions. This method + * only computes 'dg' for the reversible reactions, and the + * entries of 'dg' for the irreversible reactions are + * unaltered. This is primarily designed for use in + * calculating reveerse rate coefficients from thermochemistry + * for reversible reactions. + */ void getRevReactionDelta(int nr, const doublereal* g, doublereal* dg); - void multiplyReactants(const doublereal* c, doublereal* r); + /** + * Given an array of concentrations C, multiply the entries in array R by + * the concentration products for the reactants: + * \f[ + * R_i = R_i * \prod_k C_k^{o_{k,i}} + * \f] + * Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i. + */ + void multiplyReactants(const doublereal* C, doublereal* R); + + /** + * Given an array of concentrations C, multiply the entries in array R by + * the concentration products for the products: + * \f[ + * R_i = R_i * \prod_k C_k^{\nu^{(p)}_{k,i}} + * \f] + * Here \f$ \nu^{(p)}_{k,i} \f$ is the product stoichiometric coefficient + * of species k in reaction i. + */ void multiplyRevProducts(const doublereal* c, doublereal* r); protected: diff --git a/Cantera/src/RxnRates.h b/Cantera/src/RxnRates.h index 8ece6bdaf..b488d8663 100755 --- a/Cantera/src/RxnRates.h +++ b/Cantera/src/RxnRates.h @@ -1,7 +1,9 @@ /** * @file RxnRates.h * - * $Author$ + */ + +/* $Author$ * $Revision$ * $Date$ */ @@ -22,7 +24,7 @@ namespace Cantera { public: static int type(){ return ARRHENIUS; } Arrhenius() : m_b (0.0), m_E (0.0) {} - Arrhenius( const doublereal* c ) + Arrhenius( int csize, const doublereal* c ) : m_b (c[1]), m_E (c[2]) { m_logA = log(c[0]);} void update_C(const doublereal* c) {} @@ -46,10 +48,88 @@ namespace Cantera { return m_E; } + static bool alwaysComputeRate() { return false;} + protected: doublereal m_logA, m_b, m_E; }; + + + /** + * An Arrhenius rate with coverage-dependent terms. + */ + class SurfaceArrhenius { + + public: + static int type(){ return ARRHENIUS; } + SurfaceArrhenius() : m_b (0.0), m_E (0.0), + m_acov(0.0), m_ecov(0.0), + m_mcov(0.0), m_ncov(0), m_nmcov(0) + {} + SurfaceArrhenius( int csize, const doublereal* c ) + : m_b (c[1]), m_E (c[2]), + m_acov(0.0), m_ecov(0.0), m_mcov(0.0), m_ncov(0), m_nmcov(0) + { m_logA = log(c[0]); + if (csize >= 7) { + for (int n = 3; n < csize-3; n += 4) { + addCoverageDependence(int(c[n]), + c[n+1], c[n+2], c[n+3]); + } + } + } + + void addCoverageDependence(int k, doublereal a, + doublereal m, doublereal e) { + m_ncov++; + m_sp.push_back(k); + m_ac.push_back(a); + m_ec.push_back(e); + if (m != 0.0) { + m_msp.push_back(k); + m_mc.push_back(m); + m_nmcov++; + } + } + + void update_C(const doublereal* theta) { + m_acov = 0.0; + m_ecov = 0.0; + m_mcov = 0.0; + int n, k; + doublereal th; + for (n = 0; n < m_ncov; n++) { + k = m_sp[n]; + m_acov += m_ac[n] * theta[k]; + m_ecov += m_ec[n] * theta[k]; + } + for (n = 0; n < m_nmcov; n++) { + k = m_msp[n]; + th = fmax(theta[n], Tiny); + m_mcov += m_mc[n]*log(th); + } + } + + doublereal update(doublereal logT, doublereal recipT) const { + return m_logA + m_acov + m_b*logT + - (m_E + m_ecov)*recipT + m_mcov; + } + + doublereal activationEnergy_R() const { + return m_E + m_ecov; + } + + static bool alwaysComputeRate() { return true;} + + protected: + doublereal m_logA, m_b, m_E; + doublereal m_acov, m_ecov, m_mcov; + vector_int m_sp, m_msp; + vector_fp m_ac, m_ec, m_mc; + int m_ncov, m_nmcov; + }; + + #ifdef INCL_TST class TST { diff --git a/Cantera/src/StoichManager.h b/Cantera/src/StoichManager.h index 418109a74..7461a761f 100755 --- a/Cantera/src/StoichManager.h +++ b/Cantera/src/StoichManager.h @@ -20,8 +20,11 @@ using namespace std; namespace Cantera { - /** + /** * @defgroup Stoichiometry Stoichiometry + * + * Note: these classes are designed for internal use in class ReactionStoichManager. + * * 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 diff --git a/Cantera/src/converters/CKReader.cpp b/Cantera/src/converters/CKReader.cpp index a314f7072..8c5fda25d 100755 --- a/Cantera/src/converters/CKReader.cpp +++ b/Cantera/src/converters/CKReader.cpp @@ -216,7 +216,7 @@ bool CKReader::read(const string& inputFile, const string& thermoDatabase, } } - if (!validateSpecies(log)) { + if (validate && !validateSpecies(log)) { //Cantera::setError("read","error in species"); return false; } diff --git a/Cantera/src/converters/ck2ct.cpp b/Cantera/src/converters/ck2ct.cpp index 304b53045..2b4cf6e14 100644 --- a/Cantera/src/converters/ck2ct.cpp +++ b/Cantera/src/converters/ck2ct.cpp @@ -407,7 +407,7 @@ namespace pip { // Otherwise, just add the whole reaction, which may or may // not be reversible. - else { + else { addReaction(idktag, irxn, r.reactions[i], r.units, version); irxn++; @@ -423,7 +423,7 @@ namespace pip { const char* tr_file, const char* id_tag) { ckr::CKReader r; - r.validate = true; + r.validate = false; //int i=1; string infile = string(in_file); diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp index b24af3ed8..4bfe61d23 100755 --- a/Cantera/src/importCTML.cpp +++ b/Cantera/src/importCTML.cpp @@ -11,11 +11,6 @@ /* $Author$ * $Revision$ * $Date$ - * $Log$ - * Revision 1.14 2003-08-17 18:56:16 dggoodwin - * Added support for coverage-dependent reaction rates and sticking coefficients, - * and reactions with specified reaction order. - * */ // Copyright 2002 California Institute of Technology diff --git a/Cantera/src/importCTML.h b/Cantera/src/importCTML.h index 421fceb08..8bd23a0bc 100755 --- a/Cantera/src/importCTML.h +++ b/Cantera/src/importCTML.h @@ -14,15 +14,21 @@ #define CT_IMPORTCTML_H #include +using namespace std; -#include "Kinetics.h" -#include "transport/TransportBase.h" +#include "ThermoPhase.h" +//#include "Kinetics.h" +//#include "transport/TransportBase.h" namespace Cantera { + class Kinetics; + //class ThermoPhase; + class XML_Node; + bool isCTMLFile(string infile); - bool importPhase(XML_Node& phase, thermophase_t* th); - bool importKinetics(XML_Node& phase, vector th, + bool importPhase(XML_Node& phase, ThermoPhase* th); + bool importKinetics(XML_Node& phase, vector th, Kinetics* kin); bool installReactionArrays(XML_Node& parent, Kinetics& kin, string default_phase, bool check_for_duplicates = false); diff --git a/Cantera/src/reaction_defs.h b/Cantera/src/reaction_defs.h index 0e99dc0f8..d08c9e626 100755 --- a/Cantera/src/reaction_defs.h +++ b/Cantera/src/reaction_defs.h @@ -64,7 +64,8 @@ namespace Cantera { const int ARRHENIUS = 1; const int LANDAUTELLER = 2; const int TSTRATE = 3; - + const int SURF_ARRHENIUS = 4; + //@} /** @name Falloff Function Types