Skip to content

Commit

Permalink
Finish adding photolysis (#5)
Browse files Browse the repository at this point in the history
- Call modifyProductStoichCoeff in updateROP
- Change RateType to RateType& in MultiRate

TODO:
- [ ] Test box photolysis chemistry
  • Loading branch information
chengcli committed Apr 26, 2024
1 parent 4b813c2 commit a598366
Show file tree
Hide file tree
Showing 11 changed files with 119 additions and 3 deletions.
4 changes: 4 additions & 0 deletions include/cantera/kinetics/BulkKinetics.h
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,8 @@ class BulkKinetics : public Kinetics

void updateActinicFlux(void *rt_solver) override;

bool isPhotolysis(size_t i) const override;

//! @}

protected:
Expand Down Expand Up @@ -193,6 +195,8 @@ class BulkKinetics : public Kinetics
vector<double> m_sbuf0;
vector<double> m_state;
vector<double> m_grt; //!< Standard chemical potentials for each species

vector<size_t> m_photo_index; //!< Indices of photolysis reactions
};

}
Expand Down
6 changes: 5 additions & 1 deletion include/cantera/kinetics/Kinetics.h
Original file line number Diff line number Diff line change
Expand Up @@ -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");
}

//! @}

Expand Down
4 changes: 3 additions & 1 deletion include/cantera/kinetics/MultiRate.h
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#ifndef CT_MULTIRATE_H
#define CT_MULTIRATE_H

#include <iostream>

#include "ReactionRate.h"
#include "MultiRateBase.h"
#include "cantera/base/utilities.h"
Expand Down Expand Up @@ -203,7 +205,7 @@ class MultiRate final : public MultiRateBase
}

//! Vector of pairs of reaction rates indices and reaction rates
vector<pair<size_t, RateType>> m_rxn_rates;
vector<pair<size_t, RateType&>> m_rxn_rates;
map<size_t, size_t> m_indices; //! Mapping of indices
DataType m_shared;
};
Expand Down
2 changes: 2 additions & 0 deletions include/cantera/kinetics/MultiRateBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions include/cantera/kinetics/Photolysis.h
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
3 changes: 3 additions & 0 deletions include/cantera/kinetics/Reaction.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 5 additions & 0 deletions include/cantera/kinetics/ReactionRate.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
16 changes: 16 additions & 0 deletions src/kinetics/BulkKinetics.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
m_irrev.push_back(nReactions()-1);
}

if (r->photolysis) {
m_photo_index.push_back(nReactions()-1);
}

shared_ptr<ReactionRate> rate = r->rate();
string rtype = rate->subType();
if (rtype == "") {
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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();
}

}
2 changes: 1 addition & 1 deletion src/kinetics/Kinetics.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -693,7 +693,7 @@ shared_ptr<const Reaction> 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);

Expand Down
36 changes: 36 additions & 0 deletions src/kinetics/Reaction.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -308,6 +309,41 @@ void Reaction::setRate(shared_ptr<ReactionRate> 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<ThirdBody>("(+M)");
} else if (rate_type == "Photolysis") {
photolysis = true;
}
}
}

string Reaction::reactantString() const
Expand Down
40 changes: 40 additions & 0 deletions test/photochem/PhotochemTitan.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit a598366

Please sign in to comment.