Skip to content

Commit

Permalink
[Kinetics] Ensure rate name is known during third body set up
Browse files Browse the repository at this point in the history
  • Loading branch information
@speth
speth committed Mar 28, 2024
1 parent 5aefaef commit b0ba90a
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 62 deletions.
80 changes: 41 additions & 39 deletions src/kinetics/Reaction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,18 +64,20 @@ Reaction::Reaction(const Composition& reactants_,
m_third_body->explicit_3rd = true;
}
}
check();
}

Reaction::Reaction(const string& equation,
shared_ptr<ReactionRate> rate_,
shared_ptr<ThirdBody> tbody_)
: m_third_body(tbody_)
{
setEquation(equation);
setRate(rate_);
setEquation(equation);
if (m_third_body && m_third_body->name() != "M") {
m_third_body->explicit_3rd = true;
}
check();
}

Reaction::Reaction(const AnyMap& node, const Kinetics& kin)
Expand Down Expand Up @@ -149,10 +151,45 @@ void Reaction::check()
"Reaction orders may only be given for irreversible reactions");
}

if (!m_rate) {
return;
}

// Check reaction rate evaluator to ensure changes introduced after object
// instantiation are considered.
if (m_rate) {
m_rate->check(equation());
m_rate->check(equation());

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)");
}
}
}

Expand Down Expand Up @@ -271,39 +308,6 @@ 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)");
}
}
}

string Reaction::reactantString() const
Expand Down Expand Up @@ -354,8 +358,7 @@ string Reaction::equation() const
void Reaction::setEquation(const string& equation, const Kinetics* kin)
{
parseReactionEquation(*this, equation, input, kin);

string rate_type = input.getString("type", "");
string rate_type = (m_rate) ? m_rate->type() : input.getString("type", "");
if (ba::starts_with(rate_type, "three-body")) {
// state type when serializing
m_explicit_type = true;
Expand Down Expand Up @@ -917,7 +920,6 @@ vector<shared_ptr<Reaction>> getReactions(const AnyValue& items, Kinetics& kinet
vector<shared_ptr<Reaction>> all_reactions;
for (const auto& node : items.asVector<AnyMap>()) {
auto R = make_shared<Reaction>(node, kinetics);
R->check();
R->validate(kinetics);
if (R->valid() && R->checkSpecies(kinetics)) {
all_reactions.emplace_back(R);
Expand Down
8 changes: 2 additions & 6 deletions test/kinetics/kineticsFromScratch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,20 +486,16 @@ TEST_F(KineticsFromScratch, negative_A_error)
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
auto rate = make_shared<ArrheniusRate>(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
auto R = make_shared<Reaction>(reac, prod, rate);

ASSERT_THROW(kin.addReaction(R), CanteraError);
ASSERT_EQ((size_t) 0, kin.nReactions());
ASSERT_THROW(make_shared<Reaction>(reac, prod, rate), CanteraError);
}

TEST_F(KineticsFromScratch, allow_negative_A)
{
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
auto rate = make_shared<ArrheniusRate>(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
rate->setAllowNegativePreExponentialFactor(true);
auto R = make_shared<Reaction>(reac, prod, rate);
auto rr = std::dynamic_pointer_cast<ArrheniusRate>(R->rate());
rr->setAllowNegativePreExponentialFactor(true);

kin.addReaction(R);
ASSERT_EQ((size_t) 1, kin.nReactions());
Expand Down
4 changes: 2 additions & 2 deletions test/python/test_convert.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1465,8 +1465,8 @@ def test_stand_alone_lxcat(self):

# verify the data
assert len(rxn_list) == 2
assert rxn_list[0].equation == "O2 + e => O2(Total-Ionization)+ + e + e"
assert rxn_list[0].reaction_type == "three-body-electron-collision-plasma"
assert rxn_list[0].equation == "O2 + e => O2(Total-Ionization)+ + 2 e"
assert rxn_list[0].reaction_type == "electron-collision-plasma"
assert rxn_list[0].rate.energy_levels == approx([15., 20.])
assert rxn_list[0].rate.cross_sections == approx([0.0, 5.5e-22])

Expand Down
25 changes: 10 additions & 15 deletions test/python/test_kinetics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1283,16 +1283,14 @@ def test_arrhenius_rate(self):

def test_negative_A(self):
species = ct.Species.list_from_file("gri30.yaml")
r = ct.Reaction("NH:1, NO:1", "N2O:1, H:1",
ct.ArrheniusRate(-2.16e13, -0.23, 0))

self.assertFalse(r.rate.allow_negative_pre_exponential_factor)
rate = ct.ArrheniusRate(-2.16e13, -0.23, 0)
assert rate.allow_negative_pre_exponential_factor is False

with self.assertRaisesRegex(ct.CanteraError, 'negative pre-exponential'):
gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
species=species, reactions=[r])
r = ct.Reaction("NH:1, NO:1", "N2O:1, H:1", rate)

r.rate.allow_negative_pre_exponential_factor = True
rate.allow_negative_pre_exponential_factor = True
r = ct.Reaction("NH:1, NO:1", "N2O:1, H:1", rate)
gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
species=species, reactions=[r])

Expand Down Expand Up @@ -1501,16 +1499,13 @@ def test_BlowersMasel(self):

def test_negative_A_blowersmasel(self):
species = ct.Solution("blowers-masel.yaml").species()
r = ct.Reaction({'O':1, 'H2':1}, {'H':1, 'OH':1},
ct.BlowersMaselRate(-3.87e1, 2.7, 6260*1000*4.184, 1e9))

self.assertFalse(r.rate.allow_negative_pre_exponential_factor)

rate = ct.BlowersMaselRate(-3.87e1, 2.7, 6260*1000*4.184, 1e9)
assert rate.allow_negative_pre_exponential_factor is False
with self.assertRaisesRegex(ct.CanteraError, 'negative pre-exponential'):
gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
species=species, reactions=[r])
r = ct.Reaction({'O':1, 'H2':1}, {'H':1, 'OH':1}, rate)

r.rate.allow_negative_pre_exponential_factor = True
rate.allow_negative_pre_exponential_factor = True
r = ct.Reaction({'O':1, 'H2':1}, {'H':1, 'OH':1}, rate)
gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
species=species, reactions=[r])

Expand Down

0 comments on commit b0ba90a

Please sign in to comment.