diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 9d8af48d1e..6c62c46409 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -1625,6 +1625,9 @@ class Kinetics //! See skipUndeclaredThirdBodies() bool m_skipUndeclaredThirdBodies = false; + //! Flag indicating whether reactions include undeclared third bodies + bool m_hasUndeclaredThirdBodies = false; + //! reference to Solution std::weak_ptr m_root; }; diff --git a/interfaces/cython/cantera/kinetics.pxd b/interfaces/cython/cantera/kinetics.pxd index ca959c9451..48b2cfb23a 100644 --- a/interfaces/cython/cantera/kinetics.pxd +++ b/interfaces/cython/cantera/kinetics.pxd @@ -35,8 +35,8 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": void addThermo(shared_ptr[CxxThermoPhase]) except +translate_exception void init() except +translate_exception void skipUndeclaredThirdBodies(cbool) - void addReaction(shared_ptr[CxxReaction]) except +translate_exception - void addReaction(shared_ptr[CxxReaction], cbool) except +translate_exception + cbool addReaction(shared_ptr[CxxReaction]) except +translate_exception + cbool addReaction(shared_ptr[CxxReaction], cbool) except +translate_exception void modifyReaction(int, shared_ptr[CxxReaction]) except +translate_exception void invalidateCache() except +translate_exception void resizeReactions() diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index b4a4a4aee2..8255f502ab 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -1502,10 +1502,19 @@ cdef class Reaction: species names and the values, are the stoichiometric coefficients, for example ``{'CH4':1, 'OH':1}``, or as a composition string, for example ``'CH4:1, OH:1'``. + + .. deprecated:: 3.0 + + Setter for reactants is deprecated and will be removed after Cantera 3.0. + Use constructor instead. """ def __get__(self): return comp_map_to_dict(self.reaction.reactants) def __set__(self, reactants): + warnings.warn( + "'Reaction.reactants' setter is deprecated and will be removed after " + "Cantera 3.0.\nInstantiate using constructor instead.", + DeprecationWarning) self.reaction.reactants = comp_map(reactants) property products: @@ -1514,10 +1523,19 @@ cdef class Reaction: species names and the values, are the stoichiometric coefficients, for example ``{'CH3':1, 'H2O':1}``, or as a composition string, for example ``'CH3:1, H2O:1'``. + + .. deprecated:: 3.0 + + Setter for products is deprecated and will be removed after Cantera 3.0. + Use constructor instead. """ def __get__(self): return comp_map_to_dict(self.reaction.products) def __set__(self, products): + warnings.warn( + "'Reaction.products' setter is deprecated and will be removed after " + "Cantera 3.0.\nInstantiate using constructor instead.", + DeprecationWarning) self.reaction.products = comp_map(products) def __contains__(self, species): @@ -1655,6 +1673,18 @@ cdef class Reaction: if self.reaction.usesThirdBody(): return ThirdBody.wrap(self.reaction.thirdBody()) + @property + def third_body_name(self): + """ + Returns name of `ThirdBody` collider if `Reaction` uses a third body collider, + and ``None`` otherwise. + + .. versionadded:: 3.0 + """ + if self.reaction.usesThirdBody(): + return ThirdBody.wrap(self.reaction.thirdBody()).name + return None + property efficiencies: """ Get/Set a `dict` defining non-default third-body efficiencies for this reaction, diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index f35119164d..dfa8052ded 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -88,6 +88,8 @@ void BulkKinetics::addThirdBody(shared_ptr r) throw CanteraError("BulkKinetics::addThirdBody", "Found third-body" " efficiency for undefined species '{}' while adding reaction '{}'", name, r->equation()); + } else { + m_hasUndeclaredThirdBodies = true; } } m_multi_concm.install(nReactions() - 1, efficiencies, diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 5710d69e37..4031a705ae 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -643,6 +643,9 @@ AnyMap Kinetics::parameters() if (nReactions() == 0) { out["reactions"] = "none"; } + if (m_hasUndeclaredThirdBodies) { + out["skip-undeclared-third-bodies"] = true; + } } return out; } diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 643c83dfc1..7cb56ea1f8 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -35,21 +35,35 @@ Reaction::Reaction(const Composition& reactants_, , m_from_composition(true) , m_third_body(tbody_) { + if (reactants.count("M") || products.count("M")) { + throw CanteraError("Reaction::Reaction", + "Third body 'M' must not be included in either reactant or product maps."); + } setRate(rate_); // set flags ensuring correct serialization output - bool count = 0; + Composition third; for (const auto& [name, stoich] : reactants) { if (products.count(name)) { - count = true; - break; + third[name] = products.at(name) - stoich; } } - if (count) { - if (tbody_ && tbody_->name() != "M") { + if (tbody_) { + string name = tbody_->name(); + if (reactants.count(name) && products.count(name)) { + throw CanteraError("Reaction::Reaction", + "'{}' not acting as third body collider must not be included in both " + "reactant and product maps.", name); + } + if (name != "M") { + m_third_body->explicit_3rd = true; + } + } else if (!tbody_ && third.size() == 1) { + // implicit third body + string name = third.begin()->first; + m_third_body = make_shared(name); + if (name != "M") { m_third_body->explicit_3rd = true; - } else if (!tbody_) { - m_explicit_type = true; } } } @@ -61,7 +75,7 @@ Reaction::Reaction(const string& equation, { setEquation(equation); setRate(rate_); - if (tbody_ && tbody_->name() != "M") { + if (m_third_body && m_third_body->name() != "M") { m_third_body->explicit_3rd = true; } } @@ -822,7 +836,9 @@ bool ThirdBody::checkSpecies(const Reaction& rxn, const Kinetics& kin) const throw InputFileError("ThirdBody::checkSpecies", rxn.input, "Reaction '{}'\n" "is a three-body reaction with undeclared species: '{}'", rxn.equation(), ba::join(undeclared, "', '")); - } else if (kin.skipUndeclaredSpecies() && m_name != "M") { + } else if (kin.skipUndeclaredThirdBodies() && m_name != "M") { + // Prevent addition of reaction silently as "skip-undeclared-third-bodies" + // is set to true return false; } } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index e11304dc10..72fa122fb3 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -102,6 +102,11 @@ TEST_F(KineticsFromScratch, add_three_body_reaction1) kin.addReaction(R); check_rates(1); + + reac = parseCompString("O:2, M:1"); + prod = parseCompString("O2:1, M:1"); + ASSERT_THROW(make_shared(reac, prod, rate, tbody), CanteraError); + } TEST_F(KineticsFromScratch, add_three_body_reaction2) @@ -116,6 +121,8 @@ TEST_F(KineticsFromScratch, add_three_body_reaction2) auto tbody = make_shared(); tbody->efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4"); auto R = make_shared(equation, rate, tbody); + auto reac = R->reactants; + EXPECT_EQ(reac.count("M"), (size_t) 0); kin.addReaction(R); check_rates(1); @@ -160,6 +167,9 @@ TEST_F(KineticsFromScratch, multiple_third_bodies4) auto R = make_shared(equation, rate, tbody); EXPECT_EQ(R->thirdBody()->name(), "O2"); EXPECT_EQ(R->thirdBody()->default_efficiency, 0.); + EXPECT_EQ(R->reactants.count("H2"), (size_t) 1); + EXPECT_EQ(R->reactants.count("O2"), (size_t) 0); + EXPECT_EQ(R->reactants.count("M"), (size_t) 0); AnyMap input = R->parameters(false); EXPECT_FALSE(input.hasKey("type")); @@ -187,6 +197,9 @@ TEST_F(KineticsFromScratch, multiple_third_bodies6) auto R = make_shared(reac, prod, rate, tbody); EXPECT_EQ(R->thirdBody()->name(), "O2"); EXPECT_EQ(R->thirdBody()->default_efficiency, 0.); + EXPECT_EQ(R->reactants.count("H2"), (size_t) 1); + EXPECT_EQ(R->reactants.count("O2"), (size_t) 0); + EXPECT_EQ(R->reactants.count("M"), (size_t) 0); AnyMap input = R->parameters(false); EXPECT_FALSE(input.hasKey("type")); @@ -235,6 +248,24 @@ TEST_F(KineticsFromScratch, multiple_third_bodies8) EXPECT_EQ(input["default-efficiency"].asDouble(), 0.); } +TEST_F(KineticsFromScratch, multiple_third_bodies9) +{ + Composition reac = parseCompString("H2:1, O2:1"); + Composition prod = parseCompString("H2:1, O:2"); + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto tbody = make_shared("O2"); + auto R = make_shared(reac, prod, rate, tbody); + EXPECT_EQ(R->thirdBody()->name(), "O2"); + EXPECT_EQ(R->thirdBody()->default_efficiency, 0.); + EXPECT_EQ(R->reactants.count("H2"), (size_t) 1); + EXPECT_EQ(R->reactants.count("O2"), (size_t) 1); + EXPECT_EQ(R->reactants.count("M"), (size_t) 0); + + reac = parseCompString("H2:1, O2:1"); + prod = parseCompString("H2:1, O2:1"); + ASSERT_THROW(make_shared(reac, prod, rate, tbody), CanteraError); +} + TEST_F(KineticsFromScratch, add_two_temperature_plasma) { string equation = "O + H => O + H"; @@ -269,6 +300,35 @@ TEST_F(KineticsFromScratch, skip_undefined_third_body) ASSERT_EQ((size_t) 1, kin.nReactions()); } +TEST_F(KineticsFromScratch, skip_explicit_third_body) +{ + string equation = "2 O + CO2 <=> O2 + CO2"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto R = make_shared(equation, rate); + EXPECT_EQ(R->thirdBody()->name(), "CO2"); + + ASSERT_THROW(kin.addReaction(R), CanteraError); + kin.skipUndeclaredThirdBodies(true); + kin.addReaction(R); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch, third_body_composition) +{ + string equation = "2 O + H2O <=> O2 + H2O"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto R = make_shared(equation, rate); + EXPECT_EQ(R->thirdBody()->name(), "H2O"); + EXPECT_TRUE(R->thirdBody()->explicit_3rd); + + Composition reac = R->reactants; + EXPECT_EQ(reac.count("H2O"), (size_t) 0); + EXPECT_EQ(reac.count("M"), (size_t) 0); + Composition prod = R->products; + EXPECT_EQ(prod.count("H2O"), (size_t) 0); + EXPECT_EQ(reac.count("M"), (size_t) 0); +} + TEST_F(KineticsFromScratch, add_falloff_reaction1) { // reaction 2: diff --git a/test/python/test_kinetics.py b/test/python/test_kinetics.py index ec475ad5a0..09834ef996 100644 --- a/test/python/test_kinetics.py +++ b/test/python/test_kinetics.py @@ -304,7 +304,7 @@ def test_coverage_dependence_flags(self): surf.net_rates_of_progress_ddCi def test_electrochemistry_flags(self): - # Phases + # Phases mech = "lithium_ion_battery.yaml" anode, cathode, metal, electrolyte = ct.import_phases( mech, ["anode", "cathode", "electron", "electrolyte"]) @@ -315,6 +315,41 @@ def test_electrochemistry_flags(self): anode_int.derivative_settings = {"skip-electrochemistry": True} anode_int.net_rates_of_progress_ddCi + def test_submechanism(self): + # Simplified samples/python/kinetics/extract_submechanism.py + gri30 = ct.Solution('gri30.yaml', transport_model=None) + h2o2 = ct.Solution('h2o2.yaml', transport_model=None) + + reactions_plus = [] + reactions = [] + dest_species = set(h2o2.species_names) + colliders = dest_species.union([None, "M"]) + for R in gri30.reactions(): + if not all(S in dest_species for S in R.reactants): + continue + if not all(S in dest_species for S in R.products): + continue + reactions_plus.append(R) + if R.third_body_name not in colliders: + continue + reactions.append(R) + + gas = ct.Solution(thermo='ideal-gas', kinetics='gas', + species=h2o2.species(), reactions=reactions_plus) + # there is one third-body reaction with an undeclared third body species + gas.n_reactions < len(reactions_plus) + assert gas.n_species == len(h2o2.species_names) + + gas = ct.Solution(thermo='ideal-gas', kinetics='gas', species=h2o2.species(), + reactions=reactions) + assert gas.n_reactions == len(reactions) + assert gas.n_species == len(h2o2.species_names) + + gas.write_yaml("reduced.yaml") + restored = ct.Solution("reduced.yaml") + assert gas.species_names == restored.species_names + assert gas.reaction_equations() == restored.reaction_equations() + class KineticsRepeatability(utilities.CanteraTest): """