diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index fc9426856..b2c19e868 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -284,7 +284,7 @@ Arrhenius readArrhenius(const XML_Node& arrhenius_node) getFloat(arrhenius_node, "E", "actEnergy") / GasConstant); } -Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate_node, +Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate, const Kinetics& kin, const UnitSystem& units, int pressure_dependence=0) { @@ -314,10 +314,10 @@ Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate_node, len_dim += pressure_dependence * reaction_phase_ndim; quantity_dim -= pressure_dependence; - 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"); + auto& rate_vec = rate.asVector(); + double A = units.convert(rate_vec[0], Units(1.0, 0, len_dim, -1, 0, 0, quantity_dim)); + double b = rate_vec[1].asDouble(); + double Ta = units.convertMolarEnergy(rate_vec[2], "K"); return Arrhenius(A, b, Ta); } @@ -665,6 +665,19 @@ void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node) setupReaction(R, rxn_node); } +void setupPlogReaction(PlogReaction& R, const AnyMap& node, + const Kinetics& kin, const UnitSystem& units) +{ + setupReaction(R, node); + std::multimap rates; + for (const auto& rate : node.at("rate-constants").asVector()) { + const auto& p_rate = rate.asVector(); + rates.insert({units.convert(p_rate[0], "Pa"), + readArrhenius(R, p_rate[1], kin, units)}); + } + R.rate = Plog(rates); +} + void PlogReaction::validate() { Reaction::validate(); @@ -896,6 +909,10 @@ unique_ptr newReaction(const AnyMap& node, const Kinetics& kin, unique_ptr R(new ChemicallyActivatedReaction()); setupFalloffReaction(*R, node, kin, units); return unique_ptr(move(R)); + } else if (type == "pressure-dependent-Arrhenius") { + unique_ptr R(new PlogReaction()); + setupPlogReaction(*R, node, kin, units); + return unique_ptr(move(R)); } else { throw CanteraError("newReaction", "Unknown reaction type '{}'", type); } diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index bc0a59810..7c69d0a9e 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -117,3 +117,27 @@ TEST(Reaction, ChemicallyActivatedFromYaml) EXPECT_DOUBLE_EQ(CAR.low_rate.preExponentialFactor(), 2.82320078e2); EXPECT_EQ(CAR.falloff->nParameters(), (size_t) 0); } + +TEST(Reaction, PlogFromYaml) +{ + IdealGasMix gas("gri30.xml"); + AnyMap rxn = AnyMap::fromYamlString( + "equation: 'H + CH4 <=> H2 + CH3'\n" + "type: pressure-dependent-Arrhenius\n" + "rate-constants:\n" + "- [0.039474, [2.720000e+09 cm^3/mol/s, 1.2, 6834.0]]\n" + "- [1.0 atm, [1.260000e+20, -1.83, 15003.0]]\n" + "- [1.0 atm, [1.230000e+04, 2.68, 6335.0]]\n" + "- [1.01325 MPa, [1.680000e+16, -0.6, 14754.0]]"); + + UnitSystem U({"atm"}); + auto R = newReaction(rxn, gas, U); + auto PR = dynamic_cast(*R); + const auto& rates = PR.rate.rates(); + EXPECT_EQ(rates.size(), (size_t) 4); + EXPECT_NEAR(rates[0].first, 0.039474 * OneAtm, 1e-6); + EXPECT_NEAR(rates[2].first, OneAtm, 1e-6); + EXPECT_NEAR(rates[3].first, 10 * OneAtm, 1e-6); + EXPECT_DOUBLE_EQ(rates[0].second.preExponentialFactor(), 2.72e6); + EXPECT_DOUBLE_EQ(rates[3].second.preExponentialFactor(), 1.68e16); +}