From de80f068871a3c81ff2bcea50ca32d3683413cdc Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 5 Dec 2018 13:09:25 -0500 Subject: [PATCH] [Input] Parse YAML entries for elementary reactions --- include/cantera/base/ct_defs.h | 1 + include/cantera/kinetics/Reaction.h | 6 ++ src/kinetics/Reaction.cpp | 112 ++++++++++++++++++++++++++++ test/kinetics/kineticsFromYaml.cpp | 28 +++++++ 4 files changed, 147 insertions(+) create mode 100644 test/kinetics/kineticsFromYaml.cpp diff --git a/include/cantera/base/ct_defs.h b/include/cantera/base/ct_defs.h index 943a49ccd..ddf95f086 100644 --- a/include/cantera/base/ct_defs.h +++ b/include/cantera/base/ct_defs.h @@ -35,6 +35,7 @@ namespace Cantera using std::shared_ptr; using std::make_shared; +using std::unique_ptr; using std::isnan; // workaround for bug in libstdc++ 4.8 /*! diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 871cee1e9..cbd6529ae 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -16,6 +16,8 @@ namespace Cantera { class Kinetics; +class AnyMap; +class UnitSystem; //! Intermediate class which stores data about a reaction and its rate //! parameterization so that it can be added to a Kinetics object. @@ -257,6 +259,10 @@ public: //! Create a new Reaction object for the reaction defined in `rxn_node` shared_ptr newReaction(const XML_Node& rxn_node); +//! Create a new Reaction object using the specified parameters +unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin, + const UnitSystem& units); + //! Create Reaction objects for all `` nodes in an XML document. //! //! The `` nodes are assumed to be children of the `` diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index e121740dd..dfcb384e9 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -7,8 +7,11 @@ #include "cantera/kinetics/Reaction.h" #include "cantera/kinetics/FalloffFactory.h" +#include "cantera/kinetics/Kinetics.h" #include "cantera/base/ctml.h" #include "cantera/base/Array.h" +#include "cantera/base/AnyMap.h" +#include "cantera/base/Units.h" #include namespace Cantera @@ -277,6 +280,32 @@ Arrhenius readArrhenius(const XML_Node& arrhenius_node) getFloat(arrhenius_node, "E", "actEnergy") / GasConstant); } +Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate_node, + const Kinetics& kin, const UnitSystem& units) +{ + // Determine the units of the rate coefficient + double len_dim = - static_cast(kin.thermo(kin.reactionPhaseIndex()).nDim()); + double quantity_dim = 1.0; + for (const auto& order : R.orders) { + len_dim += order.second * kin.speciesPhase(order.first).nDim(); + quantity_dim -= order.second; + } + for (const auto& stoich : R.reactants) { + // Order for each reactant is the reactant stoichiometric coefficient, + // unless already overridden by user-specified orders + if (R.orders.find(stoich.first) == R.orders.end()) { + len_dim += stoich.second * kin.speciesPhase(stoich.first).nDim(); + quantity_dim -= stoich.second; + } + } + + const auto& rate = rate_node.asVector(); + double A = units.convert(rate[0], Units(1.0, 0, len_dim, -1, 0, 0, quantity_dim)); + double b = rate[1].asDouble(); + double Ta = units.convertMolarEnergy(rate[2], "K"); + return Arrhenius(A, b, Ta); +} + //! Parse falloff parameters, given a rateCoeff node /*! * @verbatim @@ -350,6 +379,64 @@ void setupReaction(Reaction& R, const XML_Node& rxn_node) R.reversible = (rev == "true" || rev == "yes"); } +void setupReaction(Reaction& R, const AnyMap& node) +{ + // Parse the reaction equation to determine participating species and + // stoichiometric coefficients + std::vector tokens; + tokenizeString(node.at("equation").asString(), tokens); + tokens.push_back("+"); // makes parsing last species not a special case + + size_t last_used = npos; // index of last-used token + bool reactants = true; + for (size_t i = 0; i < tokens.size(); i++) { + if (tokens[i] == "+" || tokens[i] == "<=>" || tokens[i] == "=>") { + std::string species = tokens[i-1]; + + double stoich; + if (last_used == i-2) { // Species with no stoich. coefficient + stoich = 1.0; + } else if (last_used == i-3) { // Stoich. coefficient and species + stoich = fpValueCheck(tokens[i-2]); + } else { + throw CanteraError("setupReaction", "Error parsing reaction " + "string '{}'", node.at("equation").asString()); + } + + if (reactants) { + R.reactants[species] += stoich; + } else { + R.products[species] += stoich; + } + + last_used = i; + } + + // Tokens after this point are part of the products string + if (tokens[i] == "<=>") { + R.reversible = true; + reactants = false; + } else if (tokens[i] == "=>") { + R.reversible = false; + reactants = false; + } + } + + // Non-stoichiometric reaction orders + std::map orders; + if (node.hasKey("orders")) { + for (const auto& order : node.at("orders").asMap()) { + R.orders[order.first] = order.second; + } + } + + //Flags + R.id = node.getString("id", ""); + R.duplicate = node.getBool("duplicate", false); + R.allow_negative_orders = node.getBool("negative-orders", false); + R.allow_nonreactant_orders = node.getBool("nonreactant-orders", false); +} + void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) { const XML_Node& rc_node = rxn_node.child("rateCoeff"); @@ -372,6 +459,14 @@ void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) setupReaction(R, rxn_node); } +void setupElementaryReaction(ElementaryReaction& R, const AnyMap& node, + const Kinetics& kin, const UnitSystem& units) +{ + setupReaction(R, node); + R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false); + R.rate = readArrhenius(R, node.at("rate-constant"), kin, units); +} + void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node) { readEfficiencies(R.third_body, rxn_node.child("rateCoeff")); @@ -654,6 +749,23 @@ shared_ptr newReaction(const XML_Node& rxn_node) } } +unique_ptr newReaction(const AnyMap& node, const Kinetics& kin, + const UnitSystem& units) +{ + std::string type = "elementary"; + if (node.hasKey("type")) { + type = node.at("type").asString(); + } + + if (type == "elementary") { + unique_ptr R(new ElementaryReaction()); + setupElementaryReaction(*R, node, kin, units); + return unique_ptr(move(R)); + } else { + throw CanteraError("newReaction", "Unknown reaction type '{}'", type); + } +} + std::vector > getReactions(const XML_Node& node) { std::vector > all_reactions; diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp new file mode 100644 index 000000000..0c73225aa --- /dev/null +++ b/test/kinetics/kineticsFromYaml.cpp @@ -0,0 +1,28 @@ +#include "gtest/gtest.h" +#include "cantera/base/Units.h" +#include "cantera/IdealGasMix.h" + +using namespace Cantera; + +TEST(Reaction, ElementaryFromYaml) +{ + // @TODO: Use of XML input files in these tests of the YAML format needs to + // be eliminated before we can deprecate the XML format. + IdealGasMix gas("gri30.xml"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: N + NO <=> N2 + O," + " rate-constant: [-2.70000E+13 cm^3/mol/s, 0, 355 cal/mol]," + " negative-A: true}"); + + UnitSystem U; + auto R = newReaction(rxn, gas, U); + EXPECT_EQ(R->reactants.at("NO"), 1); + EXPECT_EQ(R->products.at("N2"), 1); + EXPECT_EQ(R->reaction_type, ELEMENTARY_RXN); + + auto ER = dynamic_cast(*R); + EXPECT_DOUBLE_EQ(ER.rate.preExponentialFactor(), -2.7e10); + EXPECT_DOUBLE_EQ(ER.rate.activationEnergy_R(), 355 / GasConst_cal_mol_K); + EXPECT_TRUE(ER.allow_negative_pre_exponential_factor); + EXPECT_FALSE(ER.allow_negative_orders); +}