Skip to content

Commit

Permalink
[Kinetics] Add ReactionRate::subType to differentiate FalloffRate spe…
Browse files Browse the repository at this point in the history
…cializations
  • Loading branch information
ischoegl committed Jul 15, 2022
1 parent f6c2bdd commit fbfbc13
Show file tree
Hide file tree
Showing 11 changed files with 94 additions and 82 deletions.
13 changes: 8 additions & 5 deletions include/cantera/kinetics/Falloff.h
Expand Up @@ -170,7 +170,10 @@ class FalloffRate : public ReactionRate
}

virtual const std::string type() const {
return "Falloff";
if (m_chemicallyActivated) {
return "chemically-activated";
}
return "falloff";
}

//! Returns the number of parameters used by this parameterization. The
Expand Down Expand Up @@ -302,7 +305,7 @@ class LindemannRate final : public FalloffRate
new MultiRate<LindemannRate, FalloffData>);
}

virtual const std::string type() const {
virtual const std::string subType() const override {
return "Lindemann";
}
};
Expand Down Expand Up @@ -385,7 +388,7 @@ class TroeRate final : public FalloffRate
return 1;
}

virtual const std::string type() const {
virtual const std::string subType() const override {
return "Troe";
}

Expand Down Expand Up @@ -488,7 +491,7 @@ class SriRate final : public FalloffRate
return 2;
}

virtual const std::string type() const {
virtual const std::string subType() const override {
return "SRI";
}

Expand Down Expand Up @@ -597,7 +600,7 @@ class TsangRate final : public FalloffRate
return 1;
}

virtual const std::string type() const {
virtual const std::string subType() const override {
return "Tsang";
}

Expand Down
5 changes: 5 additions & 0 deletions include/cantera/kinetics/ReactionRate.h
Expand Up @@ -83,6 +83,11 @@ class ReactionRate
//! String identifying reaction rate specialization
virtual const std::string type() const = 0;

//! String identifying sub-type of reaction rate specialization
virtual const std::string subType() const {
return "";
}

//! Set parameters
//! @param node AnyMap object containing reaction rate specification
//! @param units unit definitions specific to rate information
Expand Down
1 change: 1 addition & 0 deletions interfaces/cython/cantera/reaction.pxd
Expand Up @@ -32,6 +32,7 @@ cdef extern from "cantera/kinetics/ReactionRate.h" namespace "Cantera":
cdef cppclass CxxReactionRate "Cantera::ReactionRate":
CxxReactionRate()
string type()
string subType()
double eval(double) except +translate_exception
double eval(double, double) except +translate_exception
double eval(double, vector[double]&) except +translate_exception
Expand Down
11 changes: 9 additions & 2 deletions interfaces/cython/cantera/reaction.pyx
Expand Up @@ -37,6 +37,11 @@ cdef class ReactionRate:
def __get__(self):
return pystr(self.rate.type())

property sub_type:
""" Get the C++ ReactionRate sub-type """
def __get__(self):
return pystr(self.rate.subType())

@staticmethod
cdef wrap(shared_ptr[CxxReactionRate] rate):
"""
Expand All @@ -53,8 +58,10 @@ cdef class ReactionRate:
# update global reaction class registry
register_subclasses(ReactionRate)

# identify class
rate_type = pystr(rate.get().type())
# identify class (either subType or type)
rate_type = pystr(rate.get().subType())
if rate_type == "":
rate_type = pystr(rate.get().type())
cls = _reaction_rate_class_registry.get(rate_type, ReactionRate)

# wrap C++ reaction rate
Expand Down
2 changes: 1 addition & 1 deletion interfaces/cython/cantera/test/test_jacobian.py
Expand Up @@ -386,7 +386,7 @@ class TestFalloff(HydrogenOxygen, utilities.CanteraTest):
# Fall-off reaction
rxn_idx = 21
equation = "2 OH (+M) <=> H2O2 (+M)"
rate_type = "Troe"
rate_type = "falloff"
rtol = 1e-4


Expand Down
6 changes: 3 additions & 3 deletions interfaces/cython/cantera/test/test_kinetics.py
Expand Up @@ -1177,7 +1177,7 @@ def test_falloff(self):
r = ct.Reaction("OH:2", "H2O2:1",
ct.TroeRate(low_rate, high_rate, [0.7346, 94, 1756, 5182]),
third_body=tb)
self.assertEqual(r.rate.type, "Troe")
self.assertEqual(r.rate.type, "falloff")

gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
species=self.species, reactions=[r])
Expand Down Expand Up @@ -1500,11 +1500,11 @@ def test_modify_falloff(self):
gas = ct.Solution('gri30.yaml', transport_model=None)
gas.TPX = 1100, 3 * ct.one_atm, 'CH4:1.0, O2:0.4, CO2:0.1, H2O:0.05'
r0 = gas.reaction(11)
self.assertEqual(r0.rate.type, "Lindemann")
self.assertEqual(r0.rate.type, "falloff")
# these two reactions happen to have the same third-body efficiencies
r1 = gas.reaction(49)
r2 = gas.reaction(53)
self.assertEqual(r2.rate.type, "Troe")
self.assertEqual(r2.rate.type, "falloff")
self.assertEqual(r1.third_body.efficiencies, r2.third_body.efficiencies)
r2.rate = r1.rate

Expand Down
11 changes: 4 additions & 7 deletions interfaces/cython/cantera/test/test_reaction.py
Expand Up @@ -409,6 +409,7 @@ def test_from_parts(self):

class FalloffRateTests(ReactionRateTests):
# test Falloff rate expressions
_type = "falloff"
_n_data = [0] # list of valid falloff coefficient array lengths

@classmethod
Expand Down Expand Up @@ -469,7 +470,6 @@ class TestLindemannRate(FalloffRateTests, utilities.CanteraTest):
# test Lindemann rate expressions

_cls = ct.LindemannRate
_type = "Lindemann"
_index = 7
_parts = {
"falloff_coeffs": [],
Expand All @@ -495,7 +495,6 @@ class TestTroeRate(FalloffRateTests, utilities.CanteraTest):
# test Troe rate expressions

_cls = ct.TroeRate
_type = "Troe"
_index = 2
_parts = {"falloff_coeffs": [0.7346, 94.0, 1756.0, 5182.0]}
_input = {
Expand Down Expand Up @@ -528,7 +527,6 @@ class TestSriRate(FalloffRateTests, utilities.CanteraTest):
# test SRI rate expressions

_cls = ct.SriRate
_type = "SRI"
_index = 8
_parts = {"falloff_coeffs": [1.1, 700.0, 1234.0, 56.0, 0.7]}
_input = {
Expand All @@ -550,7 +548,6 @@ class TestTsangRate(FalloffRateTests, utilities.CanteraTest):
# test Tsang rate expressions

_cls = ct.TsangRate
_type = "Tsang"
_index = 9
_parts = {"falloff_coeffs": [0.95, -1.0e-04]}
_input = {
Expand Down Expand Up @@ -1264,7 +1261,7 @@ class TestTroe(ReactionTests, utilities.CanteraTest):
}
_3rd_body = ct.ThirdBody("(+M)", efficiencies={"AR": 0.7, "H2": 2.0, "H2O": 6.0})
_index = 2
_rate_type = "Troe"
_rate_type = "falloff"
_yaml = """
equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 3
type: falloff
Expand Down Expand Up @@ -1302,7 +1299,7 @@ class TestLindemann(ReactionTests, utilities.CanteraTest):
}
_3rd_body = ct.ThirdBody("(+M)", efficiencies={"AR": 0.7, "H2": 2.0, "H2O": 6.0})
_index = 7
_rate_type = "Lindemann"
_rate_type = "falloff"
_yaml = """
equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 8
duplicate: true
Expand Down Expand Up @@ -1337,7 +1334,7 @@ class TestChemicallyActivated(ReactionTests, utilities.CanteraTest):
"high_P_rate_constant": {"A": 5.88E-14, "b": 6.721, "Ea": -12644997.768}
}
_index = 10
_rate_type = "Lindemann"
_rate_type = "chemically-activated"
_yaml = """
equation: H2O + OH (+M) <=> HO2 + H2 (+M) # Reaction 11
units: {length: cm, quantity: mol, activation-energy: cal/mol}
Expand Down
69 changes: 26 additions & 43 deletions interfaces/cython/cantera/yaml2ck.py
Expand Up @@ -428,69 +428,52 @@ def build_reactions_text(reactions: Iterable[ct.Reaction]):
)
)

elif reac.reaction_type.startswith("chemically-activated"):
elif reac.rate.type in ["falloff", "chemically-activated"]:
rate = reac.rate
reaction_lines.append(
arrhenius_line.format(
equation=reac.equation,
max_reaction_length=max_reaction_length,
A=rate.low_rate.pre_exponential_factor * unit_conversion_factor,
b=rate.low_rate.temperature_exponent,
E_a=rate.low_rate.activation_energy / CALORIES_CONSTANT,
)
)
unit_conversion_factor /= 1_000.0
reaction_lines.append(
high_line.format(
A=rate.high_rate.pre_exponential_factor * unit_conversion_factor,
b=rate.high_rate.temperature_exponent,
E_a=rate.high_rate.activation_energy / CALORIES_CONSTANT,
)
)
if rate.type == "Troe":
reaction_lines.append(
"TROE /"
+ " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs))
+ "/"
)
elif rate.type == "SRI":
reaction_lines.append(
"SRI /"
+ " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs))
+ "/"
)
if reac.rate.type == "falloff":
rate1 = rate.high_rate
rate2 = rate.low_rate
unit_conversion_factor2 = unit_conversion_factor * 1_000.0
other_line = low_line
else:
rate1 = rate.low_rate
rate2 = rate.high_rate
unit_conversion_factor2 = unit_conversion_factor / 1_000.0
other_line = high_line

elif reac.reaction_type.startswith("falloff"):
rate = reac.rate
reaction_lines.append(
arrhenius_line.format(
equation=reac.equation,
max_reaction_length=max_reaction_length,
A=rate.high_rate.pre_exponential_factor * unit_conversion_factor,
b=rate.high_rate.temperature_exponent,
E_a=rate.high_rate.activation_energy / CALORIES_CONSTANT,
A=rate1.pre_exponential_factor * unit_conversion_factor,
b=rate1.temperature_exponent,
E_a=rate1.activation_energy / CALORIES_CONSTANT,
)
)
unit_conversion_factor *= 1_000.0
reaction_lines.append(
low_line.format(
A=rate.low_rate.pre_exponential_factor * unit_conversion_factor,
b=rate.low_rate.temperature_exponent,
E_a=rate.low_rate.activation_energy / CALORIES_CONSTANT,
other_line.format(
A=rate2.pre_exponential_factor * unit_conversion_factor2,
b=rate2.temperature_exponent,
E_a=rate2.activation_energy / CALORIES_CONSTANT,
)
)
if rate.type == "Troe":

if reac.rate.sub_type == "Troe":
reaction_lines.append(
"TROE /"
+ " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs))
+ "/"
)
elif rate.type == "SRI":
elif reac.rate.sub_type == "SRI":
reaction_lines.append(
"SRI /"
+ " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs))
+ " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs))
+ "/"
)
elif reac.rate.sub_type != "Lindemann":
raise ValueError(
f"Unable to convert reaction type: '{reac.reaction_type}'"
)

else:
raise ValueError(f"Unknown reaction type: '{reac.reaction_type}'")
Expand Down
22 changes: 16 additions & 6 deletions src/kinetics/BulkKinetics.cpp
Expand Up @@ -136,9 +136,14 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
}

shared_ptr<ReactionRate> rate = r->rate();
std::string rtype = rate->subType();
if (rtype == "") {
rtype = rate->type();
}

// If necessary, add new MultiRate evaluator
if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) {
m_bulk_types[rate->type()] = m_bulk_rates.size();
if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
m_bulk_types[rtype] = m_bulk_rates.size();
m_bulk_rates.push_back(rate->newMultiRate());
m_bulk_rates.back()->resize(m_kk, nReactions(), nPhases());
}
Expand All @@ -148,7 +153,7 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
rate->setContext(*r, *this);

// Add reaction rate to evaluator
size_t index = m_bulk_types[rate->type()];
size_t index = m_bulk_types[rtype];
m_bulk_rates[index]->add(nReactions() - 1, *rate);

// Add reaction to third-body evaluator
Expand Down Expand Up @@ -186,14 +191,19 @@ void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
Kinetics::modifyReaction(i, rNew);

shared_ptr<ReactionRate> rate = rNew->rate();
std::string rtype = rate->subType();
if (rtype == "") {
rtype = rate->type();
}

// Ensure that MultiRate evaluator is available
if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) {
if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
throw CanteraError("BulkKinetics::modifyReaction",
"Evaluator not available for type '{}'.", rate->type());
"Evaluator not available for type '{}'.", rtype);
}

// Replace reaction rate to evaluator
size_t index = m_bulk_types[rate->type()];
size_t index = m_bulk_types[rtype];
rate->setRateIndex(i);
rate->setContext(*rNew, *this);

Expand Down
17 changes: 13 additions & 4 deletions src/kinetics/InterfaceKinetics.cpp
Expand Up @@ -428,15 +428,20 @@ bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base, bool resize)
rate->setRateIndex(nReactions() - 1);
rate->setContext(*r_base, *this);

std::string rtype = rate->subType();
if (rtype == "") {
rtype = rate->type();
}

// If necessary, add new interface MultiRate evaluator
if (m_interfaceTypes.find(rate->type()) == m_interfaceTypes.end()) {
m_interfaceTypes[rate->type()] = m_interfaceRates.size();
if (m_interfaceTypes.find(rtype) == m_interfaceTypes.end()) {
m_interfaceTypes[rtype] = m_interfaceRates.size();
m_interfaceRates.push_back(rate->newMultiRate());
m_interfaceRates.back()->resize(m_kk, nReactions(), nPhases());
}

// Add reaction rate to evaluator
size_t index = m_interfaceTypes[rate->type()];
size_t index = m_interfaceTypes[rtype];
m_interfaceRates[index]->add(nReactions() - 1, *rate);

return true;
Expand All @@ -450,7 +455,11 @@ void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
rate->setRateIndex(i);
rate->setContext(*r_base, *this);

const auto& rtype = rate->type();
std::string rtype = rate->subType();
if (rtype == "") {
rtype = rate->type();
}

// Ensure that interface MultiRate evaluator is available
if (!m_interfaceTypes.count(rtype)) {
throw CanteraError("InterfaceKinetics::modifyReaction",
Expand Down

0 comments on commit fbfbc13

Please sign in to comment.