Eliminate unnecessary specialized Reaction types

[speth]

With the introduction of the ReactionRate class, most of the behavior that is specialized for particular kinds of reactions is contained in those classes, but we still have a (more or less) 1:1 correspondence with classes derived from Reaction, even though many of them don't have any unique behavior.

Changes proposed in this pull request

• Eliminate unnecessary specializations of C++ Reaction class.

  • So far, this covers ElementaryReaction, PlogReaction, and BlowersMaselReaction, and should be easily expandable to cover CustomFunc1Reaction and TwoTempPlasmaReaction.
  • The C++ ThreeBodyReaction, FalloffReaction and ChebyshevReaction do implement some specialized behaviors, mostly around parsing and representing the reaction equation, so I think they need to stay for now.

• In the Python module, my goal is to eliminate all the specializations of Reaction, at least for homogeneous phase reactions (I think InterfaceReaction will probably need to stay)

  • This would be more or less analogous to how we handle the ThermoPhase class in Python -- with a few exceptions, there aren't wrappers for different ThermoPhase types. Right now, I'm a bit stuck based on one of the routes that's been introduced for constructing Reaction objects.

If applicable, provide an example illustrating new features this pull request is introducing

Before this PR, a PlogReaction could be created in one step as:

R = PlogReaction(
    equation="H2 + O2 <=> 2 OH",
    rate=[(1013.25, Arrhenius(1.2124e+16, -0.5779, 45491376.8)),
          (101325., Arrhenius(4.9108e+31, -4.8507, 103649395.2)),
          (1013250., Arrhenius(1.2866e+47, -9.0246, 166508556.0)),
          (10132500., Arrhenius(5.9632e+56, -11.529, 220076726.4))],
    kinetics=gas)

Or in a two steps:

R = PlogReaction({"H2": 1, "O2": 1}, {"OH": 2})
R.rate = PlogRate(
    [(1013.25, Arrhenius(1.2124e+16, -0.5779, 45491376.8)),
     (101325., Arrhenius(4.9108e+31, -4.8507, 103649395.2)),
     (1013250., Arrhenius(1.2866e+47, -9.0246, 166508556.0)),
     (10132500., Arrhenius(5.9632e+56, -11.529, 220076726.4))]

With this PR, the first example becomes:

R = Reaction(
    equation="H2 + O2 <=> 2 OH",
    rate=PlogRate([(1013.25, Arrhenius(1.2124e+16, -0.5779, 45491376.8)),
          (101325., Arrhenius(4.9108e+31, -4.8507, 103649395.2)),
          (1013250., Arrhenius(1.2866e+47, -9.0246, 166508556.0)),
          (10132500., Arrhenius(5.9632e+56, -11.529, 220076726.4))]),
    kinetics=gas)

which is pretty similar.

The second case becomes:

R = ct.Reaction({"H2": 1, "O2": 1}, {"OH": 2})
R.rate = ct.PlogRate(
    [(1013.25, ct.Arrhenius(1.2124e+16, -0.5779, 45491376.8)),
     (101325., ct.Arrhenius(4.9108e+31, -4.8507, 103649395.2)),
     (1013250., ct.Arrhenius(1.2866e+47, -9.0246, 166508556.0)),
     (10132500., ct.Arrhenius(5.9632e+56, -11.529, 220076726.4))])

This is fine in the case of PlogRate, where the base C++ reaction type is just Reaction but presents a problem for the reaction types where we still need C++ specializations of the Reaction class, as there's no information about the rate type available in the initial constructor.

What I'm wondering is whether we need this second constructor, where no rate information is provided. Would it make sense require specifying at least an empty ReactionRate object of the correct type here, e.g.

R = ct.Reaction({"H2": 1, "O2": 1}, {"OH": 2}, rate=ct.PlogRate())

I'm also wondering if there's a good way to clean up the split between specifying the reactants and products separately versus specifying the reaction equation in these constructors.

Checklist

• The pull request includes a clear description of this code change
• Commit messages have short titles and reference relevant issues
• Build passes (scons build & scons test) and unit tests address code coverage
• Style & formatting of contributed code follows contributing guidelines
• The pull request is ready for review

[ischoegl]

@speth … I ended up deleting my earlier comment as I’m curious what you’re coming up with. I believe that Reaction specializations can be removed completely in C++ relatively soon (and in Python after 2.6). From what I’m seeing over in #1181, this includes InterfaceReactions (Which means that ReactionFactory would be gone also). All we’d have to do is to track dimensionality, which is easy based on a ReactionRate when defined in the API (or phase in the context of YAML input).

[speth]

When I started this, my hope was that it would lead to enabling removal of all the Reaction specializations in C++ and along with that the ReactionFactory class, but I'm not so sure about that right now - there's a fair bit of specialization related to handling of third bodies and efficiencies for ThreeBodyReaction and FalloffReaction, and I don't know if we want to actually move that up to the base Reaction class, and I don't think it would make sense to make the ReactionRate classes know anything about the reactants / products. I'd be interested in any ideas you have that might make this possible, though.

[ischoegl]

I don't think it would make sense to make the ReactionRate classes know anything about the reactants / products.

Agreed. But for setup purposes it may be necessary to let them have a peek (without saving anything)

I'd be interested in any ideas you have that might make this possible, though.

... sure! For interface reactions, I had to do a fair amount of parsing when it comes to StickingRate (I split up SurfaceArrhenius as keeping everything together made little sense to me). The hook I came up with is

cantera/include/cantera/kinetics/ReactionRate.h

Lines 114 to 116 in d211185

	//! Provide information on associated Reaction and Kinetics
	virtual void setContext(const Reaction& rxn, const Kinetics& kin) {
	}

I am planning to use the same to eliminate ElectrochemicalReaction (and the double-parsing required for their detection).

FWIW, here's the (preliminary) implementation used for StickingRate::setContext ('lifted' from InterfaceKinetics)

cantera/src/kinetics/Coverage.cpp

Lines 176 to 245 in d211185

	void StickingRate::setContext(const Reaction& rxn, const Kinetics& kin)
	{
	setSpecies(kin.thermo().speciesNames());
	// Identify the interface phase
	size_t iInterface = npos;
	size_t min_dim = 4;
	for (size_t n = 0; n < kin.nPhases(); n++) {
	if (kin.thermo(n).nDim() < min_dim) {
	iInterface = n;
	min_dim = kin.thermo(n).nDim();
	}
	}
	std::string sticking_species = m_stickingSpecies;
	if (sticking_species == "") {
	// Identify the sticking species if not explicitly given
	bool foundStick = false;
	for (const auto& sp : rxn.reactants) {
	size_t iPhase = kin.speciesPhaseIndex(kin.kineticsSpeciesIndex(sp.first));
	if (iPhase != iInterface) {
	// Non-interface species. There should be exactly one of these
	if (foundStick) {
	throw InputFileError("StickingRate::setContext",
	rxn.input, "Multiple non-interface species ('{}' and '{}')\n"
	"found in sticking reaction: '{}'.\nSticking species "
	"must be explicitly specified.",
	sticking_species, sp.first, rxn.equation());
	}
	foundStick = true;
	sticking_species = sp.first;
	}
	}
	if (!foundStick) {
	throw InputFileError("StickingRate::setContext",
	rxn.input, "No non-interface species found "
	"in sticking reaction: '{}'", rxn.equation());
	}
	}
	m_stickingSpecies = sticking_species;
	double surface_order = 0.0;
	double multiplier = 1.0;
	// Adjust the A-factor
	for (const auto& sp : rxn.reactants) {
	size_t iPhase = kin.speciesPhaseIndex(kin.kineticsSpeciesIndex(sp.first));
	const ThermoPhase& p = kin.thermo(iPhase);
	size_t k = p.speciesIndex(sp.first);
	if (sp.first == sticking_species) {
	multiplier *= sqrt(GasConstant / (2 * Pi * p.molecularWeight(k)));
	} else {
	// Non-sticking species. Convert from coverages used in the
	// sticking probability expression to the concentration units
	// used in the mass action rate expression. For surface phases,
	// the dependence on the site density is incorporated when the
	// rate constant is evaluated, since we don't assume that the
	// site density is known at this time.
	double order = getValue(rxn.orders, sp.first, sp.second);
	const auto& surf = dynamic_cast<const SurfPhase&>(kin.thermo());
	if (&p == &surf) {
	multiplier *= pow(surf.size(k), order);
	surface_order += order;
	} else {
	multiplier *= pow(p.standardConcentration(k), -order);
	}
	}
	}
	m_surfaceOrder = surface_order;
	m_multiplier = multiplier;
	}

This function is called when the reaction is attached to a Kinetics object.

PS:

there's a fair bit of specialization related to handling of third bodies and efficiencies for ThreeBodyReaction and FalloffReaction, and I don't know if we want to actually move that up to the base Reaction class.

The situation is not too dissimilar to coverage dependence for InterfaceReaction. For that case, I ended up handling these by ReactionRate specializations.

[codecov]

Codecov Report

Merging #1183 (8a8a7bb) into main (9cee052) will increase coverage by 0.02%.
The diff coverage is 78.94%.

@@            Coverage Diff             @@
##             main    #1183      +/-   ##
==========================================
+ Coverage   65.37%   65.39%   +0.02%     
==========================================
  Files         318      318              
  Lines       46145    46109      -36     
  Branches    19615    19601      -14     
==========================================
- Hits        30166    30153      -13     
+ Misses      13475    13452      -23     
  Partials     2504     2504              

Impacted Files	Coverage Δ
include/cantera/base/AnyMap.h	84.37% <ø> (ø)
include/cantera/kinetics/Reaction.h	100.00% <ø> (ø)
include/cantera/kinetics/ReactionFactory.h	90.00% <ø> (ø)
src/thermo/ThermoFactory.cpp	71.11% <0.00%> (+0.55%)	⬆️
src/kinetics/Kinetics.cpp	72.93% <33.33%> (-0.28%)	⬇️
src/base/AnyMap.cpp	85.10% <63.63%> (+0.15%)	⬆️
src/kinetics/Reaction.cpp	80.48% <80.85%> (+1.25%)	⬆️
src/kinetics/ReactionFactory.cpp	92.07% <100.00%> (-0.04%)	⬇️
... and 1 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 9cee052...8a8a7bb. Read the comment docs.

[ischoegl]

@speth ... this looks very good as it does cut a substantial amount of clutter that accumulated as Cantera/enhancements#87 progressed: you were even able to untangle the Python API, which is something I though would have to be deferred!

I have a couple of relatively small comments. The main question I have is whether the instantiation of a ReactionRate object shouldn't be called from within Reaction.h rather than ReactionFactory, as I have some hopes that we can eliminate the latter eventually.

Regarding merges, I believe this can be safely merged around the same time as #1166, and I will build on this to rework #1181 (which would otherwise introduce a slew of new objects that would be immediately removed).

Beyond, I think it may make sense to have another look at ThirdBodyCalc as this is the only obstacle I see to removing some of the remaining reaction specializations. Of course, some changes will also be necessary for the science section on Cantera/cantera-website 😢

PS: ugh, it looks like all of my multi-line comments were converted to regular comments (I reviewed from VSCode this time).

[ischoegl]

I think that

r1 = ct.Reaction({"O":1, "H2":1}, {"H":1, "OH":1})
r1.rate = ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184)

could be combined to

r1 = ct.Reaction({"O":1, "H2":1}, {"H":1, "OH":1}, 
    ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184))

for clarity.

Same goes for instances below and similar occurrences in test_kinetics.py.

[speth]

Are you suggesting that we make rate a required argument to the Reaction constructor? I certainly think that would simplify matters.

[ischoegl]

I wasn't necessarily suggesting this, but it would probably avoid undesirable edge cases that need to be caught otherwise.

[speth]

Well, I went ahead with implementing that option, so we'll see if that ends up improving things.

[ischoegl]

I assume this is necessary as an empty type is needed to check that uninitialized rate objects work correctly?

I actually think that this PR will make uninitialized ReactionRate objects unnecessary. I.e. if no rate is attached, evaluation should fail. Likewise, generation of partial YAML input (e.g. equation without rate information) can be removed as most of the YAML content is now attached to the rate object anyways.

[speth]

Without this, the tests of the __call__ method fail for the rate types where the __call__ method takes two parameters, but the default ReactionRate object is an ArrheniusRate that only takes one.

[ischoegl]

Ah, ok. I think we can leave the handling of empty rate parameterizations for another time (but a comment should be added) … edit: actually, requiring rate as a parameter should resolve this?

[speth]

Well, it would resolve it in the sense that None would not be a viable argument to from_rate.
I put in a comment in from_rate that I hope explains this well enough. I'm looking forward to getting rid of the "legacy" rates so we can reduce the complexity of both the actual implementations and this testing infrastructure.

[ischoegl]

This leaves the object in an interesting state. Should there be any tests? I.e. what does calling the rate property getter do?

[ischoegl]

Per your comment above, requiring rate as a parameter would resolve this edge case

[speth]

Right now, this gives you a nan-initialized ArrheniusRate object. I think that ultimately stems from the fact that not specifying a reaction type implies an elementary/Arrhenius rate.

[ischoegl]

@speth ... I thought about this some more over the weekend, and ended up opening Cantera/enhancements#133.

Beyond, to unwind the remaining tasks pertaining to Cantera/enhancements#87, I believe this PR should go first, then #1184, then #1181, and finally whatever we decide for Cantera/enhancements#133. The outcome would be complete removal of Reaction specializations.

[speth]

@ischoegl - I think I've addressed all of your review comments on this.

[ischoegl]

@speth ... thank you for addressing the comments! This looks good to go - I believe #1184 should be next. As an aside, I am likewise looking forward to removing the legacy bits from the testing framework - but it does give me some peace of mind until the transition is complete.