[Input] Parse YAML entries for elementary reactions

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
 
 /*!

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<Reaction> newReaction(const XML_Node& rxn_node);
 
+//! Create a new Reaction object using the specified parameters
+unique_ptr<Reaction> newReaction(const AnyMap& rxn_node, const Kinetics& kin,
+                                 const UnitSystem& units);
+
 //! Create Reaction objects for all `<reaction>` nodes in an XML document.
 //!
 //! The `<reaction>` nodes are assumed to be children of the `<reactionData>`

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 <sstream>
 
 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<double>(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<AnyValue>();
+    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<std::string> 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<std::string, double> orders;
+    if (node.hasKey("orders")) {
+        for (const auto& order : node.at("orders").asMap<double>()) {
+            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<Reaction> newReaction(const XML_Node& rxn_node)
     }
 }
 
+unique_ptr<Reaction> 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<ElementaryReaction> R(new ElementaryReaction());
+        setupElementaryReaction(*R, node, kin, units);
+        return unique_ptr<Reaction>(move(R));
+    } else {
+        throw CanteraError("newReaction", "Unknown reaction type '{}'", type);
+    }
+}
+
 std::vector<shared_ptr<Reaction> > getReactions(const XML_Node& node)
 {
     std::vector<shared_ptr<Reaction> > all_reactions;

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<ElementaryReaction&>(*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);
+}