From fce047616eb45c6cb42c8ec82344fc8d05060f75 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 14 Apr 2021 08:27:58 -0500 Subject: [PATCH 1/7] [Kinetics] Detect three-body reactions with specified collision partners Per discussion in #967, species-specific third-body reactions should use molar concentrations for collision partner concentrations instead of activity concentrations to be consistent with three-body reactions with generic collision partners. --- include/cantera/kinetics/Reaction.h | 4 ++ .../cython/cantera/test/test_kinetics.py | 6 +- src/kinetics/Kinetics.cpp | 9 ++- src/kinetics/Reaction.cpp | 72 ++++++++++++++++--- src/kinetics/ReactionFactory.cpp | 27 +++++++ 5 files changed, 106 insertions(+), 12 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index c709b5b2b3..bbc11eb33d 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -212,6 +212,8 @@ class ThreeBodyReaction : public ElementaryReaction //! Relative efficiencies of third-body species in enhancing the reaction //! rate. ThirdBody third_body; + + bool specified = false; //!< Explicitly specified collision partner }; //! A reaction that is first-order in [M] at low pressure, like a third-body @@ -505,6 +507,8 @@ void parseReactionEquation(Reaction& R, const AnyValue& equation, const Kinetics& kin); // declarations of setup functions +void setupReaction(Reaction& R, const XML_Node& rxn_node); + void setupElementaryReaction(ElementaryReaction&, const XML_Node&); //! @internal May be changed without notice in future versions void setupElementaryReaction(ElementaryReaction&, const AnyMap&, diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 8c37c2146b..81deb2d5c3 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1221,9 +1221,11 @@ def test_modify_invalid(self): self.gas.modify_reaction(0, R2) # different reactants - R = self.gas.reaction(7) + R = self.gas.reaction(7) # implicitly defined three-body reaction + R2 = ct.ElementaryReaction(R.reactants, R.products) + R2.rate = R.rate with self.assertRaisesRegex(ct.CanteraError, 'Reactants are different'): - self.gas.modify_reaction(23, R) + self.gas.modify_reaction(23, R2) # different products R = self.gas.reaction(14) diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 97b9486f52..8dd62ad009 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -143,8 +143,13 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const continue; // No overlap in third body efficiencies } } else if (R.type() == "three-body") { - ThirdBody& tb1 = dynamic_cast(R).third_body; - ThirdBody& tb2 = dynamic_cast(other).third_body; + auto& R1 = dynamic_cast(R); + auto& R2 = dynamic_cast(other); + if (R1.specified != R2.specified) { + continue; // No overlap as one third body is explicitly specified + } + ThirdBody& tb1 = R1.third_body; + ThirdBody& tb2 = R2.third_body; bool thirdBodyOk = true; for (size_t k = 0; k < nTotalSpecies(); k++) { string s = kineticsSpeciesName(k); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index ba2899523c..c95c4a1d20 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -276,11 +276,21 @@ ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants_, } std::string ThreeBodyReaction::reactantString() const { - return ElementaryReaction::reactantString() + " + M"; + if (specified) { + return ElementaryReaction::reactantString() + " + " + + third_body.efficiencies.begin()->first; + } else { + return ElementaryReaction::reactantString() + " + M"; + } } std::string ThreeBodyReaction::productString() const { - return ElementaryReaction::productString() + " + M"; + if (specified) { + return ElementaryReaction::productString() + " + " + + third_body.efficiencies.begin()->first; + } else { + return ElementaryReaction::productString() + " + M"; + } } void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) @@ -858,6 +868,46 @@ BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate, return BlowersMasel(A, b, Ta0, w); } +bool detectEfficiencies(ThreeBodyReaction& R) +{ + for (const auto& reac : R.reactants) { + // detect explicitly specified collision partner + if (R.products.count(reac.first)) { + R.third_body.efficiencies[reac.first] = 1.; + } + } + + if (R.third_body.efficiencies.size() == 0) { + return false; + } else if (R.third_body.efficiencies.size() > 1) { + throw CanteraError("detectEfficiencies", + "Found more than one explicitly specified collision partner\n" + "in reaction '{}'.", R.equation()); + } + + R.third_body.default_efficiency = 0.; + R.specified = true; + auto sp = R.third_body.efficiencies.begin(); + + // adjust reactant coefficients + auto reac = R.reactants.find(sp->first); + if (abs(reac->second - 1) > SmallNumber) { + reac->second -= 1.; + } else { + R.reactants.erase(reac); + } + + // adjust product coefficients + auto prod = R.products.find(sp->first); + if (abs(prod->second - 1) > SmallNumber) { + prod->second -= 1.; + } else { + R.products.erase(prod); + } + + return true; +} + void setupReaction(Reaction& R, const XML_Node& rxn_node) { // Reactant and product stoichiometries @@ -1000,6 +1050,9 @@ void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node) { readEfficiencies(R.third_body, rxn_node.child("rateCoeff")); setupElementaryReaction(R, rxn_node); + if (R.third_body.efficiencies.size() == 0) { + detectEfficiencies(R); + } } void setupThreeBodyReaction(ThreeBodyReaction& R, const AnyMap& node, @@ -1007,13 +1060,16 @@ void setupThreeBodyReaction(ThreeBodyReaction& R, const AnyMap& node, { setupElementaryReaction(R, node, kin); if (R.reactants.count("M") != 1 || R.products.count("M") != 1) { - throw InputFileError("setupThreeBodyReaction", node["equation"], - "Reaction equation '{}' does not contain third body 'M'", - node["equation"].asString()); + if (!detectEfficiencies(R)) { + throw InputFileError("setupThreeBodyReaction", node["equation"], + "Reaction equation '{}' does not contain third body 'M'", + node["equation"].asString()); + } + } else { + R.reactants.erase("M"); + R.products.erase("M"); + readEfficiencies(R.third_body, node); } - R.reactants.erase("M"); - R.products.erase("M"); - readEfficiencies(R.third_body, node); } void setupFalloffReaction(FalloffReaction& R, const XML_Node& rxn_node) diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index 4ca20d400c..bdd70045a6 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -159,6 +159,17 @@ ReactionFactory::ReactionFactory() }); } +bool isThreeBody(Reaction& R) +{ + for (const auto& reac : R.reactants) { + // detect explicitly specified collision partner + if (R.products.count(reac.first)) { + return true; + } + } + return false; +} + bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) { vector_fp e_counter(kin.nPhases(), 0.0); @@ -215,6 +226,15 @@ unique_ptr newReaction(const XML_Node& rxn_node) if (type != "electrochemical") { type = R->type(); } + if (type == "elementary") { + // See if this is a three-body reaction with a specified collision partner + ElementaryReaction testReaction; + setupReaction(testReaction, rxn_node); + if (isThreeBody(testReaction)) { + type = "three-body"; + R = ReactionFactory::factory()->create(type); + } + } ReactionFactory::factory()->setup_XML(type, R, rxn_node); return unique_ptr(R); @@ -226,6 +246,13 @@ unique_ptr newReaction(const AnyMap& rxn_node, std::string type = "elementary"; if (rxn_node.hasKey("type")) { type = rxn_node["type"].asString(); + } else if (kin.thermo().nDim() == 3) { + // See if this is a three-body reaction with a specified collision partner + ElementaryReaction testReaction; + parseReactionEquation(testReaction, rxn_node["equation"], kin); + if (isThreeBody(testReaction)) { + type = "three-body"; + } } if (kin.thermo().nDim() < 3 && type == "elementary") { From 2190ac1ba10171fb50e3ae46aee3c55fb6ee8368 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 14 Apr 2021 09:23:37 -0500 Subject: [PATCH 2/7] [Kinetics] Refine criteria for the detection of third body species A reaction is identified as a three-body reaction if the following criteria are satisfied: * reaction type must not be specified * exactly one species has to occur on both reactant and product side * all stoichiometric coefficients have integer values * one side has to have exactly 3 participating species --- src/kinetics/ReactionFactory.cpp | 45 ++++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index bdd70045a6..acc2abdbeb 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -161,13 +161,40 @@ ReactionFactory::ReactionFactory() bool isThreeBody(Reaction& R) { + // detect explicitly specified collision partner + size_t found = 0; for (const auto& reac : R.reactants) { - // detect explicitly specified collision partner - if (R.products.count(reac.first)) { - return true; + auto prod = R.products.find(reac.first); + if (prod != R.products.end() && + trunc(reac.second) == reac.second && trunc(prod->second) == prod->second) { + // candidate species with integer stoichiometric coefficients on both sides + found += 1; } } - return false; + if (found != 1) { + return false; + } + + // ensure that all reactants have integer stoichiometric coefficients + size_t nreac = 0; + for (const auto& reac : R.reactants) { + if (trunc(reac.second) != reac.second) { + return false; + } + nreac += reac.second; + } + + // ensure that all products have integer stoichiometric coefficients + size_t nprod = 0; + for (const auto& prod : R.products) { + if (trunc(prod.second) != prod.second) { + return false; + } + nprod += prod.second; + } + + // either reactant or product side involves exactly three species + return (nreac == 3) || (nprod == 3); } bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) @@ -223,10 +250,8 @@ unique_ptr newReaction(const XML_Node& rxn_node) throw CanteraError("newReaction", "Unknown reaction type '" + rxn_node["type"] + "'"); } - if (type != "electrochemical") { - type = R->type(); - } - if (type == "elementary") { + if (type == "") { + // Reaction type is not specified // See if this is a three-body reaction with a specified collision partner ElementaryReaction testReaction; setupReaction(testReaction, rxn_node); @@ -235,6 +260,9 @@ unique_ptr newReaction(const XML_Node& rxn_node) R = ReactionFactory::factory()->create(type); } } + if (type != "electrochemical") { + type = R->type(); + } ReactionFactory::factory()->setup_XML(type, R, rxn_node); return unique_ptr(R); @@ -247,6 +275,7 @@ unique_ptr newReaction(const AnyMap& rxn_node, if (rxn_node.hasKey("type")) { type = rxn_node["type"].asString(); } else if (kin.thermo().nDim() == 3) { + // Reaction type is not specified // See if this is a three-body reaction with a specified collision partner ElementaryReaction testReaction; parseReactionEquation(testReaction, rxn_node["equation"], kin); From 34c882972b92745eb813675bcbe81a3ec3f3b1a2 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 14 Apr 2021 20:32:00 -0500 Subject: [PATCH 3/7] [CI] Update regression test references * test-clib: species-specific three-body reaction are reformatted * test-cxx/test-f90: reformatted species-specific three-body reactions; change of spurious net reaction rates, which should be zero * test-f77-ctlib: change of spurious net production rate for AR (gri30) --- test_problems/clib_test/output_blessed.txt | 16 +-- .../cxx_samples/cxx_demo_blessed.txt | 58 +++++----- test_problems/fortran/f77_ctlib_blessed.txt | 106 +++++++++--------- test_problems/fortran/f90_demo_blessed.txt | 62 +++++----- 4 files changed, 121 insertions(+), 121 deletions(-) diff --git a/test_problems/clib_test/output_blessed.txt b/test_problems/clib_test/output_blessed.txt index 3e9eadf687..a1ddde2d3a 100644 --- a/test_problems/clib_test/output_blessed.txt +++ b/test_problems/clib_test/output_blessed.txt @@ -107,15 +107,15 @@ CO + O2 <=> CO2 + O 1.185678e-05 CH2O + O2 <=> HCO + HO2 2.924612e-12 H + O2 + M <=> HO2 + M 4.077955e-05 - H + 2 O2 <=> HO2 + O2 9.267647e-07 - H + H2O + O2 <=> H2O + HO2 4.308072e-04 - H + N2 + O2 <=> HO2 + N2 1.003166e-04 - AR + H + O2 <=> AR + HO2 0.000000e+00 + H + O2 + O2 <=> HO2 + O2 9.267647e-07 + H + O2 + H2O <=> HO2 + H2O 4.308072e-04 + H + O2 + N2 <=> HO2 + N2 1.003166e-04 + H + O2 + AR <=> HO2 + AR 0.000000e+00 H + O2 <=> O + OH 9.566542e-02 2 H + M <=> H2 + M 3.752098e-06 - 2 H + H2 <=> 2 H2 6.185839e-08 + 2 H + H2 <=> H2 + H2 6.185839e-08 2 H + H2O <=> H2 + H2O 8.033966e-06 - CO2 + 2 H <=> CO2 + H2 1.048274e-07 + 2 H + CO2 <=> H2 + CO2 1.048274e-07 H + OH + M <=> H2O + M 3.496064e-04 H + HO2 <=> H2O + O 1.093006e-05 H + HO2 <=> H2 + O2 1.120905e-04 @@ -239,7 +239,7 @@ CH3 + CH3OH <=> CH3O + CH4 7.138914e-32 C2H4 + CH3 <=> C2H3 + CH4 3.146094e-40 C2H6 + CH3 <=> C2H5 + CH4 3.552777e-46 - H2O + HCO <=> CO + H + H2O 1.106338e-05 + HCO + H2O <=> CO + H + H2O 1.106338e-05 HCO + M <=> CO + H + M 7.088738e-06 HCO + O2 <=> CO + HO2 5.164787e-07 CH2OH + O2 <=> CH2O + HO2 1.869986e-14 @@ -376,7 +376,7 @@ CH3CHO + H => CH3 + CO + H2 1.329705e-24 CH3CHO + OH => CH3 + CO + H2O 5.882935e-24 CH3CHO + HO2 => CH3 + CO + H2O2 1.642887e-29 - CH3 + CH3CHO => CH3 + CH4 + CO 6.720371e-39 + CH3CHO + CH3 => CH4 + CO + CH3 6.720371e-39 CH2CO + H (+M) <=> CH2CHO (+M) 1.464625e-20 CH2CHO + O => CH2 + CO2 + H 1.178961e-22 CH2CHO + O2 => CH2O + CO + OH 1.668470e-25 diff --git a/test_problems/cxx_samples/cxx_demo_blessed.txt b/test_problems/cxx_samples/cxx_demo_blessed.txt index d0d5beaf7c..1b17b736f6 100644 --- a/test_problems/cxx_samples/cxx_demo_blessed.txt +++ b/test_problems/cxx_samples/cxx_demo_blessed.txt @@ -22,35 +22,35 @@ Molar Entropy: 2.4314e+05 J/kmol-K Molar cp: 32363 J/kmol-K -2 O + M <=> O2 + M 0.0056821 0.0056821 7.8063e-18 kmol/m3/s -H + O + M <=> OH + M 0.0056583 0.0056583 9.541e-18 kmol/m3/s -H2 + O <=> H + OH 134.45 134.45 5.6843e-13 kmol/m3/s -HO2 + O <=> O2 + OH 0.22106 0.22106 1.9151e-15 kmol/m3/s -H2O2 + O <=> HO2 + OH 0.014203 0.014203 -9.0206e-17 kmol/m3/s -H + O2 + M <=> HO2 + M 0.027607 0.027607 -3.1225e-17 kmol/m3/s -H + 2 O2 <=> HO2 + O2 0.014337 0.014337 -8.6736e-17 kmol/m3/s -H + H2O + O2 <=> H2O + HO2 0.64778 0.64778 2.6645e-15 kmol/m3/s -H + N2 + O2 <=> HO2 + N2 0 0 0 kmol/m3/s -AR + H + O2 <=> AR + HO2 0.069116 0.069116 2.7756e-16 kmol/m3/s -H + O2 <=> O + OH 142.55 142.55 3.1264e-13 kmol/m3/s -2 H + M <=> H2 + M 0.0015624 0.0015624 2.6021e-18 kmol/m3/s -2 H + H2 <=> 2 H2 8.8217e-05 8.8217e-05 3.2526e-19 kmol/m3/s -2 H + H2O <=> H2 + H2O 0.0054975 0.0054975 1.8215e-17 kmol/m3/s -H + OH + M <=> H2O + M 0.13038 0.13038 -4.4409e-16 kmol/m3/s -H + HO2 <=> H2O + O 0.019776 0.019776 -6.9389e-17 kmol/m3/s -H + HO2 <=> H2 + O2 0.2078 0.2078 -1.1935e-15 kmol/m3/s -H + HO2 <=> 2 OH 0.42115 0.42115 2.1649e-15 kmol/m3/s -H + H2O2 <=> H2 + HO2 0.0073139 0.0073139 -4.0766e-17 kmol/m3/s -H + H2O2 <=> H2O + OH 0.0010271 0.0010271 -5.6379e-18 kmol/m3/s -H2 + OH <=> H + H2O 229.06 229.06 -1.3927e-12 kmol/m3/s -2 OH (+M) <=> H2O2 (+M) 0.54038 0.54038 -1.3323e-15 kmol/m3/s -2 OH <=> H2O + O 375.33 375.33 -3.8085e-12 kmol/m3/s -HO2 + OH <=> H2O + O2 0.40738 0.40738 -7.7716e-16 kmol/m3/s -H2O2 + OH <=> H2O + HO2 0.00166 0.00166 -1.0192e-17 kmol/m3/s -H2O2 + OH <=> H2O + HO2 7.7286 7.7286 -4.7073e-14 kmol/m3/s -2 HO2 <=> H2O2 + O2 2.9976e-06 2.9976e-06 1.5247e-20 kmol/m3/s -2 HO2 <=> H2O2 + O2 0.00083671 0.00083671 4.12e-18 kmol/m3/s -HO2 + OH <=> H2O + O2 5.7067 5.7067 -1.1546e-14 kmol/m3/s +2 O + M <=> O2 + M 0.0056821 0.0056821 -8.6736e-18 kmol/m3/s +H + O + M <=> OH + M 0.0056583 0.0056583 8.6736e-19 kmol/m3/s +H2 + O <=> H + OH 134.45 134.45 -3.4106e-13 kmol/m3/s +HO2 + O <=> O2 + OH 0.22106 0.22106 -1.0825e-15 kmol/m3/s +H2O2 + O <=> HO2 + OH 0.014203 0.014203 -3.8164e-17 kmol/m3/s +H + O2 + M <=> HO2 + M 0.027607 0.027607 1.0755e-16 kmol/m3/s +H + O2 + O2 <=> HO2 + O2 0.014337 0.014337 5.5511e-17 kmol/m3/s +H + O2 + H2O <=> HO2 + H2O 0.64778 0.64778 2.5535e-15 kmol/m3/s +H + O2 + N2 <=> HO2 + N2 0 0 0 kmol/m3/s +H + O2 + AR <=> HO2 + AR 0.069116 0.069116 2.7756e-16 kmol/m3/s +H + O2 <=> O + OH 142.55 142.55 5.6843e-14 kmol/m3/s +2 H + M <=> H2 + M 0.0015624 0.0015624 0 kmol/m3/s +2 H + H2 <=> H2 + H2 8.8217e-05 8.8217e-05 0 kmol/m3/s +2 H + H2O <=> H2 + H2O 0.0054975 0.0054975 -8.6736e-19 kmol/m3/s +H + OH + M <=> H2O + M 0.13038 0.13038 6.9389e-16 kmol/m3/s +H + HO2 <=> H2O + O 0.019776 0.019776 7.2858e-17 kmol/m3/s +H + HO2 <=> H2 + O2 0.2078 0.2078 -3.6082e-16 kmol/m3/s +H + HO2 <=> 2 OH 0.42115 0.42115 -1.3323e-15 kmol/m3/s +H + H2O2 <=> H2 + HO2 0.0073139 0.0073139 -8.6736e-19 kmol/m3/s +H + H2O2 <=> H2O + OH 0.0010271 0.0010271 0 kmol/m3/s +H2 + OH <=> H + H2O 229.06 229.06 9.3792e-13 kmol/m3/s +2 OH (+M) <=> H2O2 (+M) 0.54038 0.54038 3.3307e-16 kmol/m3/s +2 OH <=> H2O + O 375.33 375.33 2.5011e-12 kmol/m3/s +HO2 + OH <=> H2O + O2 0.40738 0.40738 7.2164e-16 kmol/m3/s +H2O2 + OH <=> H2O + HO2 0.00166 0.00166 6.9389e-18 kmol/m3/s +H2O2 + OH <=> H2O + HO2 7.7286 7.7286 3.2863e-14 kmol/m3/s +2 HO2 <=> H2O2 + O2 2.9976e-06 2.9976e-06 -4.6587e-21 kmol/m3/s +2 HO2 <=> H2O2 + O2 0.00083671 0.00083671 -1.4095e-18 kmol/m3/s +HO2 + OH <=> H2O + O2 5.7067 5.7067 9.77e-15 kmol/m3/s Viscosity: 0.00010415 Pa-s diff --git a/test_problems/fortran/f77_ctlib_blessed.txt b/test_problems/fortran/f77_ctlib_blessed.txt index 567372a488..8d912526bd 100644 --- a/test_problems/fortran/f77_ctlib_blessed.txt +++ b/test_problems/fortran/f77_ctlib_blessed.txt @@ -1,53 +1,53 @@ - 1 1.8867924809455872E-002 107579642.8935387 - 2 1.8867924809455872E-002 -127616418.8993837 - 3 1.8867924809455872E-002 -5831548.538366192 - 4 1.8867924809455872E-002 1882808.775823964 - 5 1.8867924809455872E-002 16297577.84649889 - 6 1.8867924809455872E-002 21395925.28164294 - 7 1.8867924809455872E-002 9141390.424047336 - 8 1.8867924809455872E-002 -18212667.09474580 - 9 1.8867924809455872E-002 21655277.15429875 - 10 1.8867924809455872E-002 -3714487.701247492 - 11 1.8867924809455872E-002 -18186257.48407095 - 12 1.8867924809455872E-002 -1682714.294998080 - 13 1.8867924809455872E-002 37464941.77507096 - 14 1.8867924809455872E-002 -3392174.599151304 - 15 1.8867924809455872E-002 18021024.20126612 - 16 1.8867924809455872E-002 463920.1244349164 - 17 1.8867924809455872E-002 2806717.017385867 - 18 1.8867924809455872E-002 -1667616.834202130 - 19 1.8867924809455872E-002 1880348.120627985 - 20 1.8867924809455872E-002 -8475594.985827351 - 21 1.8867924809455872E-002 -6400868.095467533 - 22 1.8867924809455872E-002 -1406782.679495723 - 23 1.8867924809455872E-002 13553530.31813303 - 24 1.8867924809455872E-002 773602.1497574156 - 25 1.8867924809455872E-002 -1655899.720775956 - 26 1.8867924809455872E-002 -1169234.881460521 - 27 1.8867924809455872E-002 -5423450.042600359 - 28 1.8867924809455872E-002 -4477253.293424521 - 29 1.8867924809455872E-002 1678803.902351405 - 30 1.8867924809455872E-002 -7621437.918464270 - 31 1.8867924809455872E-002 2920878.667345533 - 32 1.8867924809455872E-002 1587099.067790282 - 33 1.8867924809455872E-002 -1115914.875271953 - 34 1.8867924809455872E-002 -1114194.976870384 - 35 1.8867924809455872E-002 -25381975.04578404 - 36 1.8867924809455872E-002 19123070.07861118 - 37 1.8867924809455872E-002 -5098323.231374484 - 38 1.8867924809455872E-002 -437061.2305954426 - 39 1.8867924809455872E-002 -14456790.76068574 - 40 1.8867924809455872E-002 -2054215.955030773 - 41 1.8867924809455872E-002 5671696.880337518 - 42 1.8867924809455872E-002 -877769.9776598437 - 43 1.8867924809455872E-002 -5388299.031743989 - 44 1.8867924809455872E-002 -314212.4087394696 - 45 1.8867924809455872E-002 -6148936.211230934 - 46 1.8867924809455872E-002 3674698.329455073 - 47 1.8867924809455872E-002 -2633122.304698291 - 48 1.8867924809455872E-002 31625354.14713462 - 49 1.8867924809455872E-002 -2.3592239273284576E-013 - 50 1.8867924809455872E-002 2079545.452265354 - 51 1.8867924809455872E-002 -6636261.786209078 - 52 1.8867924809455872E-002 -2239463.603775544 - 53 1.8867924809455872E-002 -528516.5887338896 + 1 1.8867924809455872E-002 107579642.89353865 + 2 1.8867924809455872E-002 -127616418.89938366 + 3 1.8867924809455872E-002 -5831548.5383661920 + 4 1.8867924809455872E-002 1882808.7758239638 + 5 1.8867924809455872E-002 16297577.846498888 + 6 1.8867924809455872E-002 21395925.281642936 + 7 1.8867924809455872E-002 9141390.4240473360 + 8 1.8867924809455872E-002 -18212667.094745796 + 9 1.8867924809455872E-002 21655277.154298745 + 10 1.8867924809455872E-002 -3714487.7012474919 + 11 1.8867924809455872E-002 -18186257.484070953 + 12 1.8867924809455872E-002 -1682714.2949980805 + 13 1.8867924809455872E-002 37464941.775070965 + 14 1.8867924809455872E-002 -3392174.5991513035 + 15 1.8867924809455872E-002 18021024.201266117 + 16 1.8867924809455872E-002 463920.12443491642 + 17 1.8867924809455872E-002 2806717.0173858674 + 18 1.8867924809455872E-002 -1667616.8342021296 + 19 1.8867924809455872E-002 1880348.1206279851 + 20 1.8867924809455872E-002 -8475594.9858273510 + 21 1.8867924809455872E-002 -6400868.0954675330 + 22 1.8867924809455872E-002 -1406782.6794957225 + 23 1.8867924809455872E-002 13553530.318133026 + 24 1.8867924809455872E-002 773602.14975741564 + 25 1.8867924809455872E-002 -1655899.7207759558 + 26 1.8867924809455872E-002 -1169234.8814605207 + 27 1.8867924809455872E-002 -5423450.0426003588 + 28 1.8867924809455872E-002 -4477253.2934245206 + 29 1.8867924809455872E-002 1678803.9023514052 + 30 1.8867924809455872E-002 -7621437.9184642704 + 31 1.8867924809455872E-002 2920878.6673455331 + 32 1.8867924809455872E-002 1587099.0677902817 + 33 1.8867924809455872E-002 -1115914.8752719532 + 34 1.8867924809455872E-002 -1114194.9768703843 + 35 1.8867924809455872E-002 -25381975.045784038 + 36 1.8867924809455872E-002 19123070.078611176 + 37 1.8867924809455872E-002 -5098323.2313744845 + 38 1.8867924809455872E-002 -437061.23059544264 + 39 1.8867924809455872E-002 -14456790.760685740 + 40 1.8867924809455872E-002 -2054215.9550307733 + 41 1.8867924809455872E-002 5671696.8803375177 + 42 1.8867924809455872E-002 -877769.97765984375 + 43 1.8867924809455872E-002 -5388299.0317439893 + 44 1.8867924809455872E-002 -314212.40873946960 + 45 1.8867924809455872E-002 -6148936.2112309337 + 46 1.8867924809455872E-002 3674698.3294550735 + 47 1.8867924809455872E-002 -2633122.3046982908 + 48 1.8867924809455872E-002 31625354.147134617 + 49 1.8867924809455872E-002 0.0000000000000000 + 50 1.8867924809455872E-002 2079545.4522653536 + 51 1.8867924809455872E-002 -6636261.7862090776 + 52 1.8867924809455872E-002 -2239463.6037755436 + 53 1.8867924809455872E-002 -528516.58873388963 diff --git a/test_problems/fortran/f90_demo_blessed.txt b/test_problems/fortran/f90_demo_blessed.txt index 2f92489cad..214a92fe4e 100644 --- a/test_problems/fortran/f90_demo_blessed.txt +++ b/test_problems/fortran/f90_demo_blessed.txt @@ -1,10 +1,10 @@ - + ******** Fortran 90 Test Program ******** Initial state properties: Temperature: 1200.0 K -Pressure: 0.10133E+06 Pa +Pressure: 0.10132E+06 Pa Density: 0.28921 kg/m3 Molar Enthalpy: 0.23514E+08 J/kmol Molar Entropy: 0.20594E+06 J/kmol-K @@ -22,35 +22,35 @@ Molar Entropy: 0.24314E+06 J/kmol-K Molar cp: 32363. J/kmol-K -2 O + M <=> O2 + M 0.56821E-02 0.56821E-02 0.78063E-17 kmol/m3/s -H + O + M <=> OH + M 0.56583E-02 0.56583E-02 0.95410E-17 kmol/m3/s -H2 + O <=> H + OH 0.13445E+03 0.13445E+03 0.56843E-12 kmol/m3/s -HO2 + O <=> O2 + OH 0.22106E+00 0.22106E+00 0.19151E-14 kmol/m3/s -H2O2 + O <=> HO2 + OH 0.14203E-01 0.14203E-01 -0.90206E-16 kmol/m3/s -H + O2 + M <=> HO2 + M 0.27607E-01 0.27607E-01 -0.31225E-16 kmol/m3/s -H + 2 O2 <=> HO2 + O2 0.14337E-01 0.14337E-01 -0.86736E-16 kmol/m3/s -H + H2O + O2 <=> H2O + HO2 0.64778E+00 0.64778E+00 0.26645E-14 kmol/m3/s -H + N2 + O2 <=> HO2 + N2 0.00000E+00 0.00000E+00 0.00000E+00 kmol/m3/s -AR + H + O2 <=> AR + HO2 0.69116E-01 0.69116E-01 0.27756E-15 kmol/m3/s -H + O2 <=> O + OH 0.14255E+03 0.14255E+03 0.31264E-12 kmol/m3/s -2 H + M <=> H2 + M 0.15624E-02 0.15624E-02 0.26021E-17 kmol/m3/s -2 H + H2 <=> 2 H2 0.88217E-04 0.88217E-04 0.32526E-18 kmol/m3/s -2 H + H2O <=> H2 + H2O 0.54975E-02 0.54975E-02 0.18215E-16 kmol/m3/s -H + OH + M <=> H2O + M 0.13038E+00 0.13038E+00 -0.44409E-15 kmol/m3/s -H + HO2 <=> H2O + O 0.19776E-01 0.19776E-01 -0.69389E-16 kmol/m3/s -H + HO2 <=> H2 + O2 0.20780E+00 0.20780E+00 -0.11935E-14 kmol/m3/s -H + HO2 <=> 2 OH 0.42115E+00 0.42115E+00 0.21649E-14 kmol/m3/s -H + H2O2 <=> H2 + HO2 0.73139E-02 0.73139E-02 -0.40766E-16 kmol/m3/s -H + H2O2 <=> H2O + OH 0.10271E-02 0.10271E-02 -0.56379E-17 kmol/m3/s -H2 + OH <=> H + H2O 0.22906E+03 0.22906E+03 -0.13927E-11 kmol/m3/s -2 OH (+M) <=> H2O2 (+M) 0.54038E+00 0.54038E+00 -0.13323E-14 kmol/m3/s -2 OH <=> H2O + O 0.37533E+03 0.37533E+03 -0.38085E-11 kmol/m3/s -HO2 + OH <=> H2O + O2 0.40738E+00 0.40738E+00 -0.77716E-15 kmol/m3/s -H2O2 + OH <=> H2O + HO2 0.16600E-02 0.16600E-02 -0.10192E-16 kmol/m3/s -H2O2 + OH <=> H2O + HO2 0.77286E+01 0.77286E+01 -0.47073E-13 kmol/m3/s -2 HO2 <=> H2O2 + O2 0.29976E-05 0.29976E-05 0.15247E-19 kmol/m3/s -2 HO2 <=> H2O2 + O2 0.83671E-03 0.83671E-03 0.41200E-17 kmol/m3/s -HO2 + OH <=> H2O + O2 0.57067E+01 0.57067E+01 -0.11546E-13 kmol/m3/s +2 O + M <=> O2 + M 0.56821E-02 0.56821E-02 -0.86736E-17 kmol/m3/s +H + O + M <=> OH + M 0.56583E-02 0.56583E-02 0.86736E-18 kmol/m3/s +H2 + O <=> H + OH 0.13445E+03 0.13445E+03 -0.34106E-12 kmol/m3/s +HO2 + O <=> O2 + OH 0.22106E+00 0.22106E+00 -0.10825E-14 kmol/m3/s +H2O2 + O <=> HO2 + OH 0.14203E-01 0.14203E-01 -0.38164E-16 kmol/m3/s +H + O2 + M <=> HO2 + M 0.27607E-01 0.27607E-01 0.10755E-15 kmol/m3/s +H + O2 + O2 <=> HO2 + O2 0.14337E-01 0.14337E-01 0.55511E-16 kmol/m3/s +H + O2 + H2O <=> HO2 + H2O 0.64778E+00 0.64778E+00 0.25535E-14 kmol/m3/s +H + O2 + N2 <=> HO2 + N2 0.00000E+00 0.00000E+00 0.00000E+00 kmol/m3/s +H + O2 + AR <=> HO2 + AR 0.69116E-01 0.69116E-01 0.27756E-15 kmol/m3/s +H + O2 <=> O + OH 0.14255E+03 0.14255E+03 0.56843E-13 kmol/m3/s +2 H + M <=> H2 + M 0.15624E-02 0.15624E-02 0.00000E+00 kmol/m3/s +2 H + H2 <=> H2 + H2 0.88217E-04 0.88217E-04 0.00000E+00 kmol/m3/s +2 H + H2O <=> H2 + H2O 0.54975E-02 0.54975E-02 -0.86736E-18 kmol/m3/s +H + OH + M <=> H2O + M 0.13038E+00 0.13038E+00 0.69389E-15 kmol/m3/s +H + HO2 <=> H2O + O 0.19776E-01 0.19776E-01 0.72858E-16 kmol/m3/s +H + HO2 <=> H2 + O2 0.20780E+00 0.20780E+00 -0.36082E-15 kmol/m3/s +H + HO2 <=> 2 OH 0.42115E+00 0.42115E+00 -0.13323E-14 kmol/m3/s +H + H2O2 <=> H2 + HO2 0.73139E-02 0.73139E-02 -0.86736E-18 kmol/m3/s +H + H2O2 <=> H2O + OH 0.10271E-02 0.10271E-02 0.00000E+00 kmol/m3/s +H2 + OH <=> H + H2O 0.22906E+03 0.22906E+03 0.93792E-12 kmol/m3/s +2 OH (+M) <=> H2O2 (+M) 0.54038E+00 0.54038E+00 0.33307E-15 kmol/m3/s +2 OH <=> H2O + O 0.37533E+03 0.37533E+03 0.25011E-11 kmol/m3/s +HO2 + OH <=> H2O + O2 0.40738E+00 0.40738E+00 0.72164E-15 kmol/m3/s +H2O2 + OH <=> H2O + HO2 0.16600E-02 0.16600E-02 0.69389E-17 kmol/m3/s +H2O2 + OH <=> H2O + HO2 0.77286E+01 0.77286E+01 0.32863E-13 kmol/m3/s +2 HO2 <=> H2O2 + O2 0.29976E-05 0.29976E-05 -0.46587E-20 kmol/m3/s +2 HO2 <=> H2O2 + O2 0.83671E-03 0.83671E-03 -0.14095E-17 kmol/m3/s +HO2 + OH <=> H2O + O2 0.57067E+01 0.57067E+01 0.97700E-14 kmol/m3/s Viscosity: 0.10415E-03 Pa-s From 3a459cb0e22ddd9450a75ab5daf546752d8552c6 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Fri, 16 Apr 2021 15:51:44 -0500 Subject: [PATCH 4/7] [Kinetics] Catch undetected third bodies in calculateRateCoeffUnits --- src/kinetics/Reaction.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index c95c4a1d20..5eccf66e6a 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -296,8 +296,17 @@ std::string ThreeBodyReaction::productString() const { void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) { ElementaryReaction::calculateRateCoeffUnits(kin); - const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex()); - rate_units *= rxn_phase.standardConcentrationUnits().pow(-1); + bool specified = false; + for (const auto& reac : reactants) { + if (reac.first != "M" && products.count(reac.first)) { + // detected specified third-body collision partner + specified = true; + } + } + if (!specified) { + const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex()); + rate_units *= rxn_phase.standardConcentrationUnits().pow(-1); + } } void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const From 33f05f85113c820cf8d1ef5664679d487d6ca7c1 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 19 Apr 2021 17:27:25 -0500 Subject: [PATCH 5/7] [CI] unit tests for detection of implicit three-body reactions --- .../cython/cantera/test/test_reaction.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index 741c5ee7e4..9377b73da6 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -4,6 +4,58 @@ from . import utilities +class TestImplicitThirdBody(utilities.CanteraTest): + + @classmethod + def setUpClass(cls): + utilities.CanteraTest.setUpClass() + cls.gas = ct.Solution('gri30.yaml') + + def test_implicit_three_body(self): + yaml1 = """ + equation: H + 2 O2 <=> HO2 + O2 + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} + """ + rxn1 = ct.Reaction.fromYaml(yaml1, self.gas) + self.assertEqual(rxn1.reaction_type, 'three-body') + + yaml2 = """ + equation: H + O2 + M <=> HO2 + M + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} + type: three-body + default-efficiency: 0 + efficiencies: {O2: 1.0} + """ + rxn2 = ct.Reaction.fromYaml(yaml2, self.gas) + self.assertEqual(rxn1.efficiencies, rxn2.efficiencies) + self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency) + + def test_non_integer_stoich(self): + yaml = """ + equation: H + 1.5 O2 <=> HO2 + O2 + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} + """ + rxn = ct.Reaction.fromYaml(yaml, self.gas) + self.assertEqual(rxn.reaction_type, 'elementary') + + def test_not_three_body(self): + yaml = """ + equation: HCNO + H <=> H + HNCO # Reaction 270 + rate-constant: {A: 2.1e+15, b: -0.69, Ea: 2850.0} + """ + rxn = ct.Reaction.fromYaml(yaml, self.gas) + self.assertEqual(rxn.reaction_type, 'elementary') + + def test_user_override(self): + yaml = """ + equation: H + 2 O2 <=> HO2 + O2 + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} + type: elementary + """ + rxn = ct.Reaction.fromYaml(yaml, self.gas) + self.assertEqual(rxn.reaction_type, 'elementary') + + class TestElementary(utilities.CanteraTest): _cls = ct.ElementaryReaction From 0ac9e8be2e6b85da2b1612f37ee2c8919826669c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 22 Apr 2021 13:49:24 -0500 Subject: [PATCH 6/7] [Kinetics] Update duplicate check for undeclared three-body reactions Also, address review comments --- include/cantera/kinetics/Reaction.h | 2 +- .../cython/cantera/test/test_kinetics.py | 6 +-- .../cython/cantera/test/test_reaction.py | 49 +++++++++++++++++-- src/kinetics/Kinetics.cpp | 9 +--- src/kinetics/Reaction.cpp | 21 +++++--- src/kinetics/ReactionFactory.cpp | 4 +- 6 files changed, 65 insertions(+), 26 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index bbc11eb33d..410839569f 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -213,7 +213,7 @@ class ThreeBodyReaction : public ElementaryReaction //! rate. ThirdBody third_body; - bool specified = false; //!< Explicitly specified collision partner + bool specified_collision_partner = false; //!< Input specifies collision partner }; //! A reaction that is first-order in [M] at low pressure, like a third-body diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 81deb2d5c3..35584571bb 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1221,11 +1221,9 @@ def test_modify_invalid(self): self.gas.modify_reaction(0, R2) # different reactants - R = self.gas.reaction(7) # implicitly defined three-body reaction - R2 = ct.ElementaryReaction(R.reactants, R.products) - R2.rate = R.rate + R = self.gas.reaction(4) with self.assertRaisesRegex(ct.CanteraError, 'Reactants are different'): - self.gas.modify_reaction(23, R2) + self.gas.modify_reaction(23, R) # different products R = self.gas.reaction(14) diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index 9377b73da6..d3f34ae0f7 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -1,4 +1,5 @@ from math import exp +from pathlib import Path import cantera as ct from . import utilities @@ -17,7 +18,9 @@ def test_implicit_three_body(self): rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ rxn1 = ct.Reaction.fromYaml(yaml1, self.gas) - self.assertEqual(rxn1.reaction_type, 'three-body') + self.assertEqual(rxn1.reaction_type, "three-body") + self.assertEqual(rxn1.default_efficiency, 0.) + self.assertEqual(rxn1.efficiencies, {"O2": 1}) yaml2 = """ equation: H + O2 + M <=> HO2 + M @@ -30,13 +33,51 @@ def test_implicit_three_body(self): self.assertEqual(rxn1.efficiencies, rxn2.efficiencies) self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency) + def test_duplicate(self): + # @todo simplify this test + # duplicates are currently only checked for import from file + gas1 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=self.gas.species(), reactions=[]) + + yaml1 = """ + equation: H + O2 + H2O <=> HO2 + H2O + rate-constant: {A: 1.126e+19, b: -0.76, Ea: 0.0} + """ + rxn1 = ct.Reaction.fromYaml(yaml1, gas1) + + yaml2 = """ + equation: H + O2 + M <=> HO2 + M + rate-constant: {A: 1.126e+19, b: -0.76, Ea: 0.0} + type: three-body + default-efficiency: 0 + efficiencies: {H2O: 1} + """ + rxn2 = ct.Reaction.fromYaml(yaml2, gas1) + + self.assertEqual(rxn1.reaction_type, rxn2.reaction_type) + self.assertEqual(rxn1.reactants, rxn2.reactants) + self.assertEqual(rxn1.products, rxn2.products) + self.assertEqual(rxn1.efficiencies, rxn2.efficiencies) + self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency) + + gas1.add_reaction(rxn1) + gas1.add_reaction(rxn2) + + fname = 'duplicate.yaml' + gas1.write_yaml(fname) + + with self.assertRaisesRegex(Exception, "Undeclared duplicate reactions"): + gas2 = ct.Solution(fname) + + Path(fname).unlink() + def test_non_integer_stoich(self): yaml = """ equation: H + 1.5 O2 <=> HO2 + O2 rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ rxn = ct.Reaction.fromYaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, 'elementary') + self.assertEqual(rxn.reaction_type, "elementary") def test_not_three_body(self): yaml = """ @@ -44,7 +85,7 @@ def test_not_three_body(self): rate-constant: {A: 2.1e+15, b: -0.69, Ea: 2850.0} """ rxn = ct.Reaction.fromYaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, 'elementary') + self.assertEqual(rxn.reaction_type, "elementary") def test_user_override(self): yaml = """ @@ -53,7 +94,7 @@ def test_user_override(self): type: elementary """ rxn = ct.Reaction.fromYaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, 'elementary') + self.assertEqual(rxn.reaction_type, "elementary") class TestElementary(utilities.CanteraTest): diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 8dd62ad009..97b9486f52 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -143,13 +143,8 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const continue; // No overlap in third body efficiencies } } else if (R.type() == "three-body") { - auto& R1 = dynamic_cast(R); - auto& R2 = dynamic_cast(other); - if (R1.specified != R2.specified) { - continue; // No overlap as one third body is explicitly specified - } - ThirdBody& tb1 = R1.third_body; - ThirdBody& tb2 = R2.third_body; + ThirdBody& tb1 = dynamic_cast(R).third_body; + ThirdBody& tb2 = dynamic_cast(other).third_body; bool thirdBodyOk = true; for (size_t k = 0; k < nTotalSpecies(); k++) { string s = kineticsSpeciesName(k); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 5eccf66e6a..c5bf263788 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -276,7 +276,7 @@ ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants_, } std::string ThreeBodyReaction::reactantString() const { - if (specified) { + if (specified_collision_partner) { return ElementaryReaction::reactantString() + " + " + third_body.efficiencies.begin()->first; } else { @@ -285,7 +285,7 @@ std::string ThreeBodyReaction::reactantString() const { } std::string ThreeBodyReaction::productString() const { - if (specified) { + if (specified_collision_partner) { return ElementaryReaction::productString() + " + " + third_body.efficiencies.begin()->first; } else { @@ -296,14 +296,19 @@ std::string ThreeBodyReaction::productString() const { void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) { ElementaryReaction::calculateRateCoeffUnits(kin); - bool specified = false; + bool specified_collision_partner = false; for (const auto& reac : reactants) { + // While this reaction was already identified as a three-body reaction in a + // pre-processing step, this method is often called before a three-body + // reaction is fully instantiated. For the determination of the correct units, + // it is necessary to check whether the reaction uses a generic 'M' or an + // explicitly specified collision partner that may not have been deleted yet. if (reac.first != "M" && products.count(reac.first)) { // detected specified third-body collision partner - specified = true; + specified_collision_partner = true; } } - if (!specified) { + if (!specified_collision_partner) { const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex()); rate_units *= rxn_phase.standardConcentrationUnits().pow(-1); } @@ -895,12 +900,12 @@ bool detectEfficiencies(ThreeBodyReaction& R) } R.third_body.default_efficiency = 0.; - R.specified = true; + R.specified_collision_partner = true; auto sp = R.third_body.efficiencies.begin(); // adjust reactant coefficients auto reac = R.reactants.find(sp->first); - if (abs(reac->second - 1) > SmallNumber) { + if (trunc(reac->second) != 1) { reac->second -= 1.; } else { R.reactants.erase(reac); @@ -908,7 +913,7 @@ bool detectEfficiencies(ThreeBodyReaction& R) // adjust product coefficients auto prod = R.products.find(sp->first); - if (abs(prod->second - 1) > SmallNumber) { + if (trunc(prod->second) != 1) { prod->second -= 1.; } else { R.products.erase(prod); diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index acc2abdbeb..b0b7834ac6 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -159,7 +159,7 @@ ReactionFactory::ReactionFactory() }); } -bool isThreeBody(Reaction& R) +bool isThreeBody(const Reaction& R) { // detect explicitly specified collision partner size_t found = 0; @@ -250,7 +250,7 @@ unique_ptr newReaction(const XML_Node& rxn_node) throw CanteraError("newReaction", "Unknown reaction type '" + rxn_node["type"] + "'"); } - if (type == "") { + if (type.empty()) { // Reaction type is not specified // See if this is a three-body reaction with a specified collision partner ElementaryReaction testReaction; From f61f85c0351ce00daed7f5bd9798ee4d9c996fbf Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 22 Apr 2021 15:37:46 -0500 Subject: [PATCH 7/7] [Kinetics] Update serialization of three-body reactions Three-body reactions with explictly declared collision partners do not need to declare additional fields, i.e. 'type', 'default-efficiency', and 'efficiencies'. --- .../cython/cantera/test/test_reaction.py | 18 +++++++++++++++--- src/kinetics/Reaction.cpp | 18 ++++++++++-------- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index d3f34ae0f7..47e54d85d7 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -10,7 +10,7 @@ class TestImplicitThirdBody(utilities.CanteraTest): @classmethod def setUpClass(cls): utilities.CanteraTest.setUpClass() - cls.gas = ct.Solution('gri30.yaml') + cls.gas = ct.Solution("gri30.yaml") def test_implicit_three_body(self): yaml1 = """ @@ -36,7 +36,7 @@ def test_implicit_three_body(self): def test_duplicate(self): # @todo simplify this test # duplicates are currently only checked for import from file - gas1 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + gas1 = ct.Solution(thermo="IdealGas", kinetics="GasKinetics", species=self.gas.species(), reactions=[]) yaml1 = """ @@ -63,7 +63,7 @@ def test_duplicate(self): gas1.add_reaction(rxn1) gas1.add_reaction(rxn2) - fname = 'duplicate.yaml' + fname = "duplicate.yaml" gas1.write_yaml(fname) with self.assertRaisesRegex(Exception, "Undeclared duplicate reactions"): @@ -71,6 +71,18 @@ def test_duplicate(self): Path(fname).unlink() + def test_short_serialization(self): + yaml = """ + equation: H + O2 + H2O <=> HO2 + H2O + rate-constant: {A: 1.126e+19, b: -0.76, Ea: 0.0} + """ + rxn = ct.Reaction.fromYaml(yaml, self.gas) + input_data = rxn.input_data + + self.assertNotIn("type", input_data) + self.assertNotIn("default-efficiency", input_data) + self.assertNotIn("efficiencies", input_data) + def test_non_integer_stoich(self): yaml = """ equation: H + 1.5 O2 <=> HO2 + O2 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index c5bf263788..165ab3af3d 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -296,7 +296,7 @@ std::string ThreeBodyReaction::productString() const { void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) { ElementaryReaction::calculateRateCoeffUnits(kin); - bool specified_collision_partner = false; + bool specified_collision_partner_ = false; for (const auto& reac : reactants) { // While this reaction was already identified as a three-body reaction in a // pre-processing step, this method is often called before a three-body @@ -305,10 +305,10 @@ void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) // explicitly specified collision partner that may not have been deleted yet. if (reac.first != "M" && products.count(reac.first)) { // detected specified third-body collision partner - specified_collision_partner = true; + specified_collision_partner_ = true; } } - if (!specified_collision_partner) { + if (!specified_collision_partner_) { const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex()); rate_units *= rxn_phase.standardConcentrationUnits().pow(-1); } @@ -317,11 +317,13 @@ void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const { ElementaryReaction::getParameters(reactionNode); - reactionNode["type"] = "three-body"; - reactionNode["efficiencies"] = third_body.efficiencies; - reactionNode["efficiencies"].setFlowStyle(); - if (third_body.default_efficiency != 1.0) { - reactionNode["default-efficiency"] = third_body.default_efficiency; + if (!specified_collision_partner) { + reactionNode["type"] = "three-body"; + reactionNode["efficiencies"] = third_body.efficiencies; + reactionNode["efficiencies"].setFlowStyle(); + if (third_body.default_efficiency != 1.0) { + reactionNode["default-efficiency"] = third_body.default_efficiency; + } } }