From a5983666b0eaeacddd02986071b03824f9c4fded Mon Sep 17 00:00:00 2001 From: Cheng Li <69489965+chengcli@users.noreply.github.com> Date: Thu, 4 Jan 2024 09:17:10 -0500 Subject: [PATCH] Finish adding photolysis (#5) - Call modifyProductStoichCoeff in updateROP - Change RateType to RateType& in MultiRate TODO: - [ ] Test box photolysis chemistry --- include/cantera/kinetics/BulkKinetics.h | 4 +++ include/cantera/kinetics/Kinetics.h | 6 +++- include/cantera/kinetics/MultiRate.h | 4 ++- include/cantera/kinetics/MultiRateBase.h | 2 ++ include/cantera/kinetics/Photolysis.h | 4 +++ include/cantera/kinetics/Reaction.h | 3 ++ include/cantera/kinetics/ReactionRate.h | 5 +++ src/kinetics/BulkKinetics.cpp | 16 ++++++++++ src/kinetics/Kinetics.cpp | 2 +- src/kinetics/Reaction.cpp | 36 +++++++++++++++++++++ test/photochem/PhotochemTitan.cpp | 40 ++++++++++++++++++++++++ 11 files changed, 119 insertions(+), 3 deletions(-) diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index c67abb3f17..44224f7495 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -98,6 +98,8 @@ class BulkKinetics : public Kinetics void updateActinicFlux(void *rt_solver) override; + bool isPhotolysis(size_t i) const override; + //! @} protected: @@ -193,6 +195,8 @@ class BulkKinetics : public Kinetics vector m_sbuf0; vector m_state; vector m_grt; //!< Standard chemical potentials for each species + + vector m_photo_index; //!< Indices of photolysis reactions }; } diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 2d66d595b0..759b1b080d 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -1450,7 +1450,11 @@ class Kinetics * The set of products must not be changed. * Only the stoichiometric coefficient is changed. */ - virtual void modifyProductStoichiometry(size_t i, Composition const& comp); + virtual void modifyProductStoichCoeff(size_t i, Composition const& comp); + + virtual bool isPhotolysis(size_t i) const { + throw NotImplementedError("Kinetics::isPhotolysis"); + } //! @} diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h index 48aeb1c9cd..6e8ddc315d 100644 --- a/include/cantera/kinetics/MultiRate.h +++ b/include/cantera/kinetics/MultiRate.h @@ -8,6 +8,8 @@ #ifndef CT_MULTIRATE_H #define CT_MULTIRATE_H +#include + #include "ReactionRate.h" #include "MultiRateBase.h" #include "cantera/base/utilities.h" @@ -203,7 +205,7 @@ class MultiRate final : public MultiRateBase } //! Vector of pairs of reaction rates indices and reaction rates - vector> m_rxn_rates; + vector> m_rxn_rates; map m_indices; //! Mapping of indices DataType m_shared; }; diff --git a/include/cantera/kinetics/MultiRateBase.h b/include/cantera/kinetics/MultiRateBase.h index 33d4182944..8a67c91699 100644 --- a/include/cantera/kinetics/MultiRateBase.h +++ b/include/cantera/kinetics/MultiRateBase.h @@ -32,6 +32,8 @@ class MultiRateBase //! Identifier of reaction rate type virtual string type() = 0; + virtual void peek() const {}; + //! Add reaction rate object to the evaluator //! @param rxn_index index of reaction //! @param rate reaction rate object diff --git a/include/cantera/kinetics/Photolysis.h b/include/cantera/kinetics/Photolysis.h index 3d70edd4f3..e9aed1d12b 100644 --- a/include/cantera/kinetics/Photolysis.h +++ b/include/cantera/kinetics/Photolysis.h @@ -124,6 +124,10 @@ class PhotolysisRate : public PhotolysisBase { return "Photolysis"; } + Composition const& photoProducts() const override { + return m_net_products; + } + double evalFromStruct(PhotolysisData const& data) { double wmin = m_temp_wave_grid[m_ntemp]; double wmax = m_temp_wave_grid.back(); diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index bd148dea15..7acbbf4767 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -162,6 +162,9 @@ class Reaction return bool(m_third_body); } + //! True if the current reaction is photolysis. False otherwise + bool photolysis = false; + protected: //! Store the parameters of a Reaction needed to reconstruct an identical //! object using the newReaction(AnyMap&, Kinetics&) function. Does not diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index b2a9bc976d..517da388e6 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -216,6 +216,11 @@ class ReactionRate m_composition_dependent_rate = comp_dep; } + virtual Composition const& photoProducts() const { + throw NotImplementedError("ReactionRate::photoProducts", + "Not implemented by '{}' object.", type()); + } + protected: //! Get parameters //! @param node AnyMap containing rate information diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 9be7631c5b..444da4cbe3 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -39,6 +39,10 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) m_irrev.push_back(nReactions()-1); } + if (r->photolysis) { + m_photo_index.push_back(nReactions()-1); + } + shared_ptr rate = r->rate(); string rtype = rate->subType(); if (rtype == "") { @@ -536,6 +540,14 @@ void BulkKinetics::updateROP() AssertFinite(m_ropr[i], "BulkKinetics::updateROP", "m_ropr[{}] is not finite.", i); } + + for (auto i : m_photo_index) { + // modify reaction if there is more than one photolysis product + if (m_reactions[i]->products.size() > 1) { + modifyProductStoichCoeff(i, m_reactions[i]->rate()->photoProducts()); + } + } + m_ROP_ok = true; } @@ -707,4 +719,8 @@ void __attribute__((weak)) BulkKinetics::updateActinicFlux(void *rt_solver) m_hasNewActinicFlux = true; } +bool BulkKinetics::isPhotolysis(size_t i) const { + return std::find(m_photo_index.begin(), m_photo_index.end(), i) < m_photo_index.end(); +} + } diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index c6218cdf08..d02b6e5d5a 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -693,7 +693,7 @@ shared_ptr Kinetics::reaction(size_t i) const return m_reactions[i]; } -void Kinetics::modifyProductStoichiometry(size_t irxn, Composition const& comp) +void Kinetics::modifyProductStoichCoeff(size_t irxn, Composition const& comp) { checkReactionIndex(irxn); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index ad13285183..6fe5ad3c0f 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -74,6 +74,7 @@ Reaction::Reaction(const string& equation, { setRate(rate_); setEquation(equation); + if (m_third_body && m_third_body->name() != "M") { m_third_body->explicit_3rd = true; } @@ -308,6 +309,41 @@ void Reaction::setRate(shared_ptr rate) "Reaction rate for reaction '{}' must not be empty.", equation()); } m_rate = rate; + + string rate_type = m_rate->type(); + if (m_third_body) { + if (rate_type == "falloff" || rate_type == "chemically-activated") { + if (m_third_body->mass_action && !m_from_composition) { + throw InputFileError("Reaction::setRate", input, + "Third-body collider does not use '(+{})' notation.", + m_third_body->name()); + } + m_third_body->mass_action = false; + } else if (rate_type == "Chebyshev") { + if (m_third_body->name() == "M") { + warn_deprecated("Chebyshev reaction equation", input, "Specifying 'M' " + "in the reaction equation for Chebyshev reactions is deprecated."); + m_third_body.reset(); + } + } else if (rate_type == "pressure-dependent-Arrhenius") { + if (m_third_body->name() == "M") { + throw InputFileError("Reaction::setRate", input, + "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.", + m_third_body->name()); + } + } + } else { + if (rate_type == "falloff" || rate_type == "chemically-activated") { + if (!m_from_composition) { + throw InputFileError("Reaction::setRate", input, + "Reaction equation for falloff reaction '{}'\n does not " + "contain valid pressure-dependent third body", equation()); + } + m_third_body = make_shared("(+M)"); + } else if (rate_type == "Photolysis") { + photolysis = true; + } + } } string Reaction::reactantString() const diff --git a/test/photochem/PhotochemTitan.cpp b/test/photochem/PhotochemTitan.cpp index 0a70858104..a57c76e306 100644 --- a/test/photochem/PhotochemTitan.cpp +++ b/test/photochem/PhotochemTitan.cpp @@ -53,6 +53,46 @@ TEST_F(PhotochemTitan, check_fwd_rate_constants) { ASSERT_NEAR(kfwd[0], 3.06820e-14, 1.0e-18); ASSERT_NEAR(kfwd[1], 3.2e-16, 1.0e-18); + + int iCH4 = kin->kineticsSpeciesIndex("CH4"); + int iCH3 = kin->kineticsSpeciesIndex("CH3"); + int i1CH2 = kin->kineticsSpeciesIndex("(1)CH2"); + int i3CH2 = kin->kineticsSpeciesIndex("(3)CH2"); + int iCH = kin->kineticsSpeciesIndex("CH"); + int iH2 = kin->kineticsSpeciesIndex("H2"); + int iH = kin->kineticsSpeciesIndex("H"); + + ASSERT_EQ(iCH4, 0); + ASSERT_EQ(iCH3, 1); + ASSERT_EQ(i1CH2, 2); + ASSERT_EQ(i3CH2, 3); + ASSERT_EQ(iCH, 4); + ASSERT_EQ(iH2, 5); + ASSERT_EQ(iH, 6); + + double kCH4 = kin->productStoichCoeff(iCH4, 0); + ASSERT_NEAR(kCH4, 0.635657, 1.0e-4); + + double kCH3 = kin->productStoichCoeff(iCH3, 0); + ASSERT_NEAR(kCH3, 0.142168, 1.0e-4); + + double k1CH2 = kin->productStoichCoeff(i1CH2, 0); + ASSERT_NEAR(k1CH2, 0.0978033, 1.0e-4); + + double k3CH2 = kin->productStoichCoeff(i3CH2, 0); + ASSERT_NEAR(k3CH2, 0.0377844, 1.0e-4); + + double kCH = kin->productStoichCoeff(iCH, 0); + ASSERT_NEAR(kCH, 0.0865869, 1.0e-4); + + double kH2 = kin->productStoichCoeff(iH2, 0); + ASSERT_NEAR(kH2, 0.18439, 1.0e-4); + + double kH = kin->productStoichCoeff(iH, 0); + ASSERT_NEAR(kH, 0.304324, 1.0e-4); + + ASSERT_NEAR(kCH4 + kCH3 + k1CH2 + k3CH2 + kCH, 1.0, 1.0e-14); + ASSERT_NEAR(4 * kCH4 + 3 * kCH3 + 2 * k1CH2 + 2 * k3CH2 + kCH + 2 * kH2 + kH, 4.0, 1.0e-14); } } // namespace Cantera