diff --git a/doc/sphinx/cti/reactions.rst b/doc/sphinx/cti/reactions.rst
index f713f3fb9..6bd361f77 100644
--- a/doc/sphinx/cti/reactions.rst
+++ b/doc/sphinx/cti/reactions.rst
@@ -92,6 +92,8 @@ Note that the ID string is only used when selectively importing reactions. If
all reactions in the local file or in an external one are imported into a phase
or interface, then the reaction ``ID`` field is not used.
+.. _sec-reaction-options:
+
Options
-------
@@ -136,6 +138,12 @@ should be handled.
positive, then negative *A* parameters are acceptable, as long as the
``'negative_A'`` option is specified.
+``negative_orders``
+ Reaction orders are normally required to be non-negative, since negative
+ orders are non-physical and undefined at zero concentration. Cantera allows
+ negative orders for a global reaction only if the ``negative_orders``
+ override option is specified for the reaction.
+
Reactions with Pressure-Independent Rate
========================================
diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h
index 38ba218e3..b9510c88b 100644
--- a/include/cantera/kinetics/Reaction.h
+++ b/include/cantera/kinetics/Reaction.h
@@ -28,8 +28,9 @@ public:
virtual std::string productString() const;
std::string equation() const;
- //! Ensure that the rate constant for this reaction is valid.
- virtual void validateRateConstant() {}
+ //! Ensure that the rate constant and other parameters for this reaction are
+ //valid.
+ virtual void validate();
//! Type of the reaction. The valid types are listed in the file,
//! reaction_defs.h, with constants ending in `RXN`.
@@ -55,6 +56,13 @@ public:
//! True if the current reaction is marked as duplicate
bool duplicate;
+
+ //! True if reaction orders can be specified for non-reactant species.
+ //Default is `false`.
+ bool allow_nonreactant_orders;
+
+ //! True if negative reaction orders are allowed. Default is `false`.
+ bool allow_negative_orders;
};
@@ -66,7 +74,7 @@ public:
ElementaryReaction();
ElementaryReaction(const Composition& reactants, const Composition products,
const Arrhenius& rate);
- virtual void validateRateConstant();
+ virtual void validate();
Arrhenius rate;
bool allow_negative_pre_exponential_factor;
@@ -115,7 +123,7 @@ public:
const vector_fp& falloff_params);
virtual std::string reactantString() const;
virtual std::string productString() const;
- virtual void validateRateConstant();
+ virtual void validate();
Arrhenius low_rate;
Arrhenius high_rate;
@@ -148,7 +156,7 @@ public:
PlogReaction();
PlogReaction(const Composition& reactants, const Composition& products,
const Plog& rate);
- virtual void validateRateConstant();
+ virtual void validate();
Plog rate;
};
diff --git a/interfaces/cython/cantera/ctml_writer.py b/interfaces/cython/cantera/ctml_writer.py
index f0bd8bf84..47ee358b0 100644
--- a/interfaces/cython/cantera/ctml_writer.py
+++ b/interfaces/cython/cantera/ctml_writer.py
@@ -1128,8 +1128,9 @@ class reaction(object):
An optional identification string. If omitted, it defaults to a
four-digit numeric string beginning with 0001 for the first
reaction in the file.
- :param options:
- Processing options, as described in :ref:`sec-phase-options`.
+ :param options: Processing options, as described in
+ :ref:`sec-reaction-options`. May be one or more (as a list) of the
+ following: 'skip', 'duplicate', 'negative_A', 'negative_orders'.
"""
self._id = id
self._e = equation
@@ -1218,11 +1219,12 @@ class reaction(object):
else:
r['reversible'] = 'no'
- for s in self._options:
- if s == 'duplicate':
- r['duplicate'] = 'yes'
- elif s == 'negative_A':
- r['negative_A'] = 'yes'
+ if 'duplicate' in self._options:
+ r['duplicate'] = 'yes'
+ if 'negative_A' in self._options:
+ r['negative_A'] = 'yes'
+ if 'negative_orders' in self._options:
+ r['negative_orders'] = 'yes'
ee = self._e.replace('<','[').replace('>',']')
r.addChild('equation',ee)
@@ -1337,8 +1339,8 @@ class three_body_reaction(reaction):
An optional identification string. If omitted, it defaults to a
four-digit numeric string beginning with 0001 for the first
reaction in the file.
- :param options:
- Processing options, as described in :ref:`sec-phase-options`.
+ :param options: Processing options, as described in
+ :ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, kf, id, '', options)
self._type = 'threeBody'
@@ -1426,7 +1428,7 @@ class falloff_reaction(pdep_reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
- Processing options, as described in :ref:`sec-phase-options`.
+ Processing options, as described in :ref:`sec-reaction-options`.
"""
kf2 = (kf, kf0)
reaction.__init__(self, equation, kf2, id, '', options)
@@ -1469,7 +1471,7 @@ class chemically_activated_reaction(pdep_reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
- Processing options, as described in :ref:`sec-phase-options`.
+ Processing options, as described in :ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, (kLow, kHigh), id, '', options)
self._type = 'chemAct'
@@ -1599,7 +1601,7 @@ class surface_reaction(reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
- Processing options, as described in :ref:`sec-phase-options`.
+ Processing options, as described in :ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, kf, id, order, options)
self._type = 'surface'
diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp
index 339853b8c..96a68b429 100644
--- a/src/kinetics/Kinetics.cpp
+++ b/src/kinetics/Kinetics.cpp
@@ -655,7 +655,7 @@ void Kinetics::addReaction(ReactionData& r) {
void Kinetics::addReaction(shared_ptr r)
{
- r->validateRateConstant();
+ r->validate();
// If reaction orders are specified, then this reaction does not follow
// mass-action kinetics, and is not an elementary reaction. So check that it
diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp
index fd1acc62f..e98470a42 100644
--- a/src/kinetics/Reaction.cpp
+++ b/src/kinetics/Reaction.cpp
@@ -16,6 +16,8 @@ Reaction::Reaction(int type)
: reaction_type(type)
, reversible(true)
, duplicate(false)
+ , allow_nonreactant_orders(false)
+ , allow_negative_orders(false)
{
}
@@ -26,9 +28,36 @@ Reaction::Reaction(int type, const Composition& reactants_,
, products(products_)
, reversible(true)
, duplicate(false)
+ , allow_nonreactant_orders(false)
+ , allow_negative_orders(false)
{
}
+void Reaction::validate()
+{
+ if (!allow_nonreactant_orders) {
+ for (Composition::iterator iter = orders.begin();
+ iter != orders.end();
+ ++iter) {
+ if (reactants.find(iter->first) == reactants.end()) {
+ throw CanteraError("Reaction::validate", "Reaction order "
+ "specified for non-reactant species '" + iter->first + "'");
+ }
+ }
+ }
+
+ if (!allow_negative_orders) {
+ for (Composition::iterator iter = orders.begin();
+ iter != orders.end();
+ ++iter) {
+ if (iter->second < 0.0) {
+ throw CanteraError("Reaction::validate", "Negative reaction "
+ "order specified for species '" + iter->first + "'");
+ }
+ }
+ }
+}
+
std::string Reaction::reactantString() const
{
std::ostringstream result;
@@ -87,11 +116,12 @@ ElementaryReaction::ElementaryReaction()
{
}
-void ElementaryReaction::validateRateConstant()
+void ElementaryReaction::validate()
{
+ Reaction::validate();
if (!allow_negative_pre_exponential_factor &&
rate.preExponentialFactor() < 0) {
- throw CanteraError("ElementaryReaction::validateRateConstant",
+ throw CanteraError("ElementaryReaction::validate",
"Undeclared negative pre-exponential factor found in reaction '"
+ equation() + "'");
}
@@ -165,10 +195,11 @@ std::string FalloffReaction::productString() const {
}
}
-void FalloffReaction::validateRateConstant() {
+void FalloffReaction::validate() {
+ Reaction::validate();
if (low_rate.preExponentialFactor() < 0 ||
high_rate.preExponentialFactor() < 0) {
- throw CanteraError("FalloffReaction::validateRateConstant", "Negative "
+ throw CanteraError("FalloffReaction::validate", "Negative "
"pre-exponential factor found for reaction '" + equation() + "'");
}
}
@@ -342,6 +373,9 @@ void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node)
if (rxn_node["negative_A"] == "yes") {
R.allow_negative_pre_exponential_factor = true;
}
+ if (rxn_node["negative_orders"] == "yes") {
+ R.allow_negative_orders = true;
+ }
setupReaction(R, rxn_node);
}
@@ -421,8 +455,9 @@ void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node)
setupReaction(R, rxn_node);
}
-void PlogReaction::validateRateConstant()
+void PlogReaction::validate()
{
+ Reaction::validate();
rate.validate(equation());
}
@@ -508,6 +543,7 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R,
R.orders.clear();
// Reaction orders based on species stoichiometric coefficients
+ R.allow_nonreactant_orders = true;
for (Composition::const_iterator iter = R.reactants.begin();
iter != R.reactants.end();
++iter) {
@@ -524,6 +560,7 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R,
if (rxn_node.hasChild("reactionOrderFormulation")) {
Composition initial_orders = R.orders;
R.orders.clear();
+ R.allow_nonreactant_orders = true;
const XML_Node& rof_node = rxn_node.child("reactionOrderFormulation");
if (lowercase(rof_node["model"]) == "reactantorders") {
R.orders = initial_orders;
diff --git a/test/data/frac.cti b/test/data/frac.cti
index a338a42b3..e78717965 100644
--- a/test/data/frac.cti
+++ b/test/data/frac.cti
@@ -97,3 +97,7 @@ reaction( "H2O => 1.4 H + 0.6 OH + 0.2 O2", [1.0e13, 0.0, 0.0])
# coefficients.
reaction( "0.7 H2 + 0.6 OH + 0.2 O2 => H2O", [1.0e13, 0.0, 0.0],
order = "H2:0.8 OH:2 O2:1")
+
+# A reaction with negative reaction orders
+reaction( "H2 + 0.5 O2 => H2O", [1.0e9, 0.0, 0.0],
+ order = "H2:1.0 O2:-0.25", options='negative_orders')
diff --git a/test/data/frac.xml b/test/data/frac.xml
index 8b09e3c3c..53a0b40b5 100644
--- a/test/data/frac.xml
+++ b/test/data/frac.xml
@@ -153,5 +153,21 @@
H2:0.69999999999999996 O2:0.20000000000000001 OH:0.59999999999999998
H2O:1.0
+
+
+ H2 + 0.5 O2 =] H2O
+ 1.0
+ -0.25
+
+
+ 3.981072E+00
+ 0.0
+ 0.000000
+
+
+ H2:1.0 O2:0.5
+ H2O:1.0
+
+
diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp
index ceacc2078..bcf1f789e 100644
--- a/test/kinetics/kineticsFromScratch.cpp
+++ b/test/kinetics/kineticsFromScratch.cpp
@@ -261,6 +261,59 @@ TEST_F(KineticsFromScratch, invalid_reversible_with_orders)
ASSERT_EQ(0, kin.nReactions());
}
+TEST_F(KineticsFromScratch, negative_order_override)
+{
+ Composition reac = parseCompString("O:1 H2:1");
+ Composition prod = parseCompString("H:1 OH:1");
+ Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
+ shared_ptr R(new ElementaryReaction(reac, prod, rate));
+ R->reversible = false;
+ R->allow_negative_orders = true;
+ R->orders["H2"] = - 0.5;
+
+ kin.addReaction(R);
+ ASSERT_EQ((size_t) 1, kin.nReactions());
+}
+
+TEST_F(KineticsFromScratch, invalid_negative_orders)
+{
+ Composition reac = parseCompString("O:1 H2:1");
+ Composition prod = parseCompString("H:1 OH:1");
+ Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
+ shared_ptr R(new ElementaryReaction(reac, prod, rate));
+ R->reversible = false;
+ R->orders["H2"] = - 0.5;
+
+ ASSERT_THROW(kin.addReaction(R), CanteraError);
+ ASSERT_EQ(0, kin.nReactions());
+}
+
+TEST_F(KineticsFromScratch, nonreactant_order_override)
+{
+ Composition reac = parseCompString("O:1 H2:1");
+ Composition prod = parseCompString("H:1 OH:1");
+ Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
+ shared_ptr R(new ElementaryReaction(reac, prod, rate));
+ R->reversible = false;
+ R->allow_nonreactant_orders = true;
+ R->orders["OH"] = 0.5;
+
+ kin.addReaction(R);
+ ASSERT_EQ((size_t) 1, kin.nReactions());
+}
+
+TEST_F(KineticsFromScratch, invalid_nonreactant_order)
+{
+ Composition reac = parseCompString("O:1 H2:1");
+ Composition prod = parseCompString("H:1 OH:1");
+ Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
+ shared_ptr R(new ElementaryReaction(reac, prod, rate));
+ R->reversible = false;
+ R->orders["OH"] = 0.5;
+
+ ASSERT_THROW(kin.addReaction(R), CanteraError);
+ ASSERT_EQ(0, kin.nReactions());
+}
class InterfaceKineticsFromScratch : public testing::Test
{
diff --git a/test/kinetics/rates.cpp b/test/kinetics/rates.cpp
index 7256d4967..384d8a839 100644
--- a/test/kinetics/rates.cpp
+++ b/test/kinetics/rates.cpp
@@ -116,10 +116,10 @@ TEST_F(FracCoeffTest, CreationDestructionRates)
EXPECT_DOUBLE_EQ(0.6*ropf[0], cdot[kOH]);
EXPECT_DOUBLE_EQ(0.2*ropf[0], cdot[kO2]);
- EXPECT_DOUBLE_EQ(0.7*ropf[1], ddot[kH2]);
+ EXPECT_DOUBLE_EQ(0.7*ropf[1]+ropf[2], ddot[kH2]);
EXPECT_DOUBLE_EQ(0.6*ropf[1], ddot[kOH]);
- EXPECT_DOUBLE_EQ(0.2*ropf[1], ddot[kO2]);
- EXPECT_DOUBLE_EQ(ropf[1], cdot[kH2O]);
+ EXPECT_DOUBLE_EQ(0.2*ropf[1]+0.5*ropf[2], ddot[kO2]);
+ EXPECT_DOUBLE_EQ(ropf[1]+ropf[2], cdot[kH2O]);
EXPECT_DOUBLE_EQ(0.0, cdot[therm.speciesIndex("O")]);
EXPECT_DOUBLE_EQ(0.0, ddot[therm.speciesIndex("O")]);