From b4c07978e6baf1d0ea382e9467646cb0068c5599 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 30 Mar 2012 23:47:31 +0000 Subject: [PATCH] Implemented parsing of P-log and Chebyshev reactions from CTML --- include/cantera/kinetics/GasKinetics.h | 2 + include/cantera/kinetics/reaction_defs.h | 28 +++- src/kinetics/GasKinetics.cpp | 29 +++- src/kinetics/ReactionData.h | 15 ++ src/kinetics/importKinetics.cpp | 179 ++++++++++++----------- 5 files changed, 156 insertions(+), 97 deletions(-) diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index dab35e75a..c397e55e7 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -452,6 +452,8 @@ private: void addElementaryReaction(const ReactionData& r); void addThreeBodyReaction(const ReactionData& r); void addFalloffReaction(const ReactionData& r); + void addPlogReaction(const ReactionData& r); + void addChebyshevReaction(const ReactionData& r); void installReagents(const ReactionData& r); diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index c03f60c45..df902ae0a 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -37,6 +37,20 @@ const int THREE_BODY_RXN = 2; */ const int FALLOFF_RXN = 4; +/** + * A pressure-dependent rate expression consisting of several Arrhenius rate + * expressions evaluated at different pressures. The final rate is calculated + * by logarithmically interpolating between the two rates that bracket the + * current pressure. + */ +const int PLOG_RXN = 5; + +/** + * A general pressure-dependent reaction where k(T,P) is defined in terms of + * a bivariate Chebyshev polynomial. + */ +const int CHEBYSHEV_RXN = 6; + /** * A chemical activation reaction. For these reactions, the rate falls * off as the pressure increases, due to collisional stabilization of @@ -75,12 +89,14 @@ const int GLOBAL_RXN = 30; */ //@{ -const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1; -const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2; -const int TSTRATE_REACTION_RATECOEFF_TYPE = 3; -const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4; -const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5; -const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6; +const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1; +const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2; +const int TSTRATE_REACTION_RATECOEFF_TYPE = 3; +const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4; +const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5; +const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6; +const int PLOG_REACTION_RATECOEFF_TYPE = 7; +const int CHEBYSHEV_REACTION_RATECOEFF_TYPE = 8; //@} diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 562d71851..73ac5ca9c 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -657,13 +657,24 @@ getRevRateConstants(doublereal* krev, bool doIrreversible) void GasKinetics:: addReaction(const ReactionData& r) { - - if (r.reactionType == ELEMENTARY_RXN) { + switch (r.reactionType) { + case ELEMENTARY_RXN: addElementaryReaction(r); - } else if (r.reactionType == THREE_BODY_RXN) { + break; + case THREE_BODY_RXN: addThreeBodyReaction(r); - } else if (r.reactionType == FALLOFF_RXN) { + break; + case FALLOFF_RXN: addFalloffReaction(r); + break; + case PLOG_RXN: + addPlogReaction(r); + break; + case CHEBYSHEV_RXN: + addChebyshevReaction(r); + break; + default: + throw CanteraError("GasKinetics::addReaction", "Invalid reaction type specified"); } // operations common to all reaction types @@ -761,6 +772,16 @@ addThreeBodyReaction(const ReactionData& r) } //==================================================================================================================== +void GasKinetics::addPlogReaction(const ReactionData& r) +{ + // @todo: Not yet implemented +} + +void GasKinetics::addChebyshevReaction(const ReactionData& r) +{ + // @todo: Not yet implemented +} + void GasKinetics::installReagents(const ReactionData& r) { diff --git a/src/kinetics/ReactionData.h b/src/kinetics/ReactionData.h index c84e02c91..567625568 100644 --- a/src/kinetics/ReactionData.h +++ b/src/kinetics/ReactionData.h @@ -9,6 +9,7 @@ #define CT_REACTION_DATA_H #include "cantera/kinetics/reaction_defs.h" +#include "cantera/kinetics/RxnRates.h" namespace Cantera { @@ -71,6 +72,20 @@ public: bool global; bool isReversibleWithFrac; doublereal beta; // for electrochemical reactions + + //! Arrhenius parameters for P-log reactions. + //! The keys are the pressures corresponding to each Arrhenius expression. + std::map plogParameters; + + double chebTmin; //!< Minimum temperature for Chebyshev fit + double chebTmax; //!< Maximum temperature for Chebyshev fit + double chebPmin; //!< Minimum pressure for Chebyshev fit + double chebPmax; //!< Maximum pressure for Chebyshev fit + size_t chebDegreeT; //!< Degree of Chebyshev fit in T + size_t chebDegreeP; //!< Degree of Chebyshev fit in P + + //! Chebyshev coefficients. length chebDegreeT * chebDegreeP + vector_fp chebCoeffs; }; } diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index dff5369b1..1243e9732 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -68,19 +68,6 @@ public: bool validate_rxn) ; }; - - -/* - * First we define a couple of typedefs that will - * be used throught this file - */ -//! typedef for a pointer to an XML_Node -typedef const vector nodeset_t; -//! typedef for an XML_Node -typedef XML_Node node_t; - - - /* * Check a reaction to see if the elements balance. */ @@ -395,7 +382,7 @@ static void getStick(const XML_Node& node, Kinetics& kin, E /= GasConstant; } -static void getCoverageDependence(const node_t& node, +static void getCoverageDependence(const XML_Node& node, thermo_t& surfphase, ReactionData& rdata) { vector cov; @@ -428,7 +415,7 @@ static void getCoverageDependence(const node_t& node, 0.5 73.2 5000. 9999. @endverbatim */ -static void getFalloff(const node_t& f, ReactionData& rdata) +static void getFalloff(const XML_Node& f, ReactionData& rdata) { string type = f["type"]; vector p; @@ -489,7 +476,7 @@ static void getFalloff(const node_t& f, ReactionData& rdata) * reaction mechanism is homogeneous, so that all species belong * to phase(0) of 'kin'. */ -static void getEfficiencies(const node_t& eff, Kinetics& kin, ReactionData& rdata) +static void getEfficiencies(const XML_Node& eff, Kinetics& kin, ReactionData& rdata) { // set the default collision efficiency @@ -514,89 +501,106 @@ static void getEfficiencies(const node_t& eff, Kinetics& kin, ReactionData& rdat * * @param kf Reference to the XML Node named rateCoeff */ -void getRateCoefficient(const node_t& kf, Kinetics& kin, +void getRateCoefficient(const XML_Node& kf, Kinetics& kin, ReactionData& rdata, int negA) { - string type = kf.attrib("type"); - if (type == "") { - type = "Arrhenius"; - rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE; - } - if (type == "ExchangeCurrentDensity") { - rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; - } else if (type == "Arrhenius") { + if (rdata.reactionType == PLOG_RXN) { + rdata.rateCoeffType = PLOG_REACTION_RATECOEFF_TYPE; + for (size_t m = 0; m < kf.nChildren(); m++) { + const XML_Node& node = kf.child(m); + double A = getFloat(node, "A", "toSI"); + double b = getFloat(node, "b"); + double E = getFloat(node, "E", "actEnergy") / GasConstant; + double p = getFloat(node, "P", "toSI"); + rdata.plogParameters[p] = Arrhenius(A, b, E); + } + + } else if (rdata.reactionType == CHEBYSHEV_RXN) { + rdata.rateCoeffType = CHEBYSHEV_REACTION_RATECOEFF_TYPE; + rdata.chebTmin = getFloat(kf, "Tmin", "toSI"); + rdata.chebTmax = getFloat(kf, "Tmax", "toSI"); + rdata.chebPmin = getFloat(kf, "Pmin", "toSI"); + rdata.chebPmax = getFloat(kf, "Pmax", "toSI"); + const XML_Node& coeffs = kf.child("floatArray"); + rdata.chebDegreeP = atoi(coeffs["degreeP"].c_str()); + rdata.chebDegreeT = atoi(coeffs["degreeT"].c_str()); + getFloatArray(kf, rdata.chebCoeffs, false); } else { - throw CanteraError("getRateCoefficient", - "Unknown type: " + type); - } - nodeset_t& kf_children = kf.children(); - vector_fp clow(3,0.0), chigh(3,0.0); - for (size_t m = 0; m < kf.nChildren(); m++) { - const node_t& c = *kf_children[m]; - string nm = c.name(); - int highlow=0; + string type = kf.attrib("type"); + if (type == "") { + type = "Arrhenius"; + rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE; + } + if (type == "ExchangeCurrentDensity") { + rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; + } else if (type == "Arrhenius") { - if (nm == "Arrhenius") { - vector_fp coeff(3); - if (c["type"] == "stick") { - getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]); - chigh = coeff; - } else { - getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); - if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN - || rdata.reactionType == ELEMENTARY_RXN) { + } else { + throw CanteraError("getRateCoefficient", "Unknown type: " + type); + } + + vector_fp clow(3,0.0), chigh(3,0.0); + for (size_t m = 0; m < kf.nChildren(); m++) { + const XML_Node& c = kf.child(m); + string nm = c.name(); + int highlow=0; + + if (nm == "Arrhenius") { + vector_fp coeff(3); + if (c["type"] == "stick") { + getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]); chigh = coeff; } else { - clow = coeff; + getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); + if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN + || rdata.reactionType == ELEMENTARY_RXN) { + chigh = coeff; + } else { + clow = coeff; + } + } + if (rdata.reactionType == SURFACE_RXN) { + getCoverageDependence(c, + kin.thermo(kin.surfacePhaseIndex()), rdata); } - } - if (rdata.reactionType == SURFACE_RXN) { - getCoverageDependence(c, - kin.thermo(kin.surfacePhaseIndex()), rdata); - } - if (coeff[0] <= 0.0 && negA == 0) { - throw CanteraError("getRateCoefficient", - "negative or zero A coefficient for reaction "+int2str(rdata.number)); + if (coeff[0] <= 0.0 && negA == 0) { + throw CanteraError("getRateCoefficient", + "negative or zero A coefficient for reaction "+int2str(rdata.number)); + } + } else if (nm == "Arrhenius_ExchangeCurrentDensity") { + vector_fp coeff(3); + getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); + chigh = coeff; + rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; + } else if (nm == "falloff") { + getFalloff(c, rdata); + } else if (nm == "efficiencies") { + getEfficiencies(c, kin, rdata); + } else if (nm == "electrochem") { + rdata.beta = fpValue(c["beta"]); } - } else if (nm == "Arrhenius_ExchangeCurrentDensity") { - vector_fp coeff(3); - getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); - chigh = coeff; - rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; - } else if (nm == "falloff") { - getFalloff(c, rdata); - } else if (nm == "efficiencies") { - getEfficiencies(c, kin, rdata); - } else if (nm == "electrochem") { - rdata.beta = fpValue(c["beta"]); + } + /* + * Store the coefficients in the ReactionData object for return + * from this function. + */ + if (rdata.reactionType == CHEMACT_RXN) { + rdata.rateCoeffParameters = clow; + } else { + rdata.rateCoeffParameters = chigh; + } + + if (rdata.reactionType == FALLOFF_RXN) { + rdata.auxRateCoeffParameters = clow; + } else if (rdata.reactionType == CHEMACT_RXN) { + rdata.auxRateCoeffParameters = chigh; } } - /* - * Store the coefficients in the ReactionData object for return - * from this function. - */ - if (rdata.reactionType == CHEMACT_RXN) { - rdata.rateCoeffParameters = clow; - } else { - rdata.rateCoeffParameters = chigh; - } - - if (rdata.reactionType == FALLOFF_RXN) { - rdata.auxRateCoeffParameters = clow; - } else if (rdata.reactionType == CHEMACT_RXN) { - rdata.auxRateCoeffParameters = chigh; - } - - - } - - - /* * This function returns true if two reactions are duplicates of * one another, and false otherwise. The input arguments are two @@ -842,6 +846,10 @@ bool rxninfo::installReaction(int i, const XML_Node& r, Kinetics* k, rdata.falloffType = SIMPLE_FALLOFF; } else if (typ == "threeBody") { rdata.reactionType = THREE_BODY_RXN; + } else if (typ == "plog") { + rdata.reactionType = PLOG_RXN; + } else if (typ == "chebyshev") { + rdata.reactionType = CHEBYSHEV_RXN; } else if (typ == "surface") { rdata.reactionType = SURFACE_RXN; } else if (typ == "edge") { @@ -1208,7 +1216,4 @@ bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm, } } - - } -