Eliminate unnecessary specialized Reaction types by speth · Pull Request #1183 · Cantera/cantera

speth · 2022-01-24T22:33:13Z

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 · 2022-01-25T01:20:13Z

@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 · 2022-01-25T02:08:30Z

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 · 2022-01-25T02:26:38Z

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 · 2022-01-26T05:25:56Z

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).

test/data/kineticsfromscratch.yaml

src/kinetics/Reaction.cpp

src/kinetics/ReactionFactory.cpp

ischoegl · 2022-01-27T17:40:15Z

interfaces/cython/cantera/examples/kinetics/blowers_masel.py


 #Create an elementary reaction O+H2<=>H+OH
-r1 = ct.ElementaryReaction({'O':1, 'H2':1}, {'H':1, 'OH':1})
+r1 = ct.Reaction({"O":1, "H2":1}, {"H":1, "OH":1})


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.

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

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

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

interfaces/cython/cantera/test/test_reaction.py

ischoegl · 2022-01-27T17:54:11Z

interfaces/cython/cantera/test/test_reaction.py

        if self._rate_obj is None:
            return
-        rxn = self.from_rate(None)
+        rxn = self.from_rate(self._rate_cls() if self._rate_cls else None)


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.

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.

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?

~~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.

interfaces/cython/cantera/reaction.pyx

ischoegl · 2022-01-27T18:03:32Z

interfaces/cython/cantera/reaction.pyx

+            return
+
+        if reactants and products and not equation:
+            equation = self.equation


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

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

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.

interfaces/cython/cantera/reaction.pyx

ischoegl · 2022-01-31T14:55:41Z

@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.

An elementary reaction is just a Reaction with an ArrheniusRate as its rate coefficient.

This is useful for cases where the deprecation refers to something that appears in the input file, rather than being a deprecated method.

While these reactions are pressure dependent, that dependence is fully incorporated into the rate coefficient and the concentration of third bodies is never directly used in calculating the reaction rate.

This maintains the standing behavior of treating different reaction types as non-duplicates after making reactions with different rate types use the base Reaction class.

speth · 2022-02-04T17:54:35Z

@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.

ischoegl mentioned this pull request Jan 25, 2022

Electrochemical Reactions missing from Science Section Cantera/cantera-website#172

Closed

speth force-pushed the simplify-reactions branch from a4a2373 to 97c20b6 Compare January 26, 2022 05:03

speth force-pushed the simplify-reactions branch 3 times, most recently from e1a4d74 to 0ec03e3 Compare January 27, 2022 04:25

speth marked this pull request as ready for review January 27, 2022 04:27

ischoegl self-assigned this Jan 27, 2022

ischoegl requested changes Jan 27, 2022

View reviewed changes

This was referenced Jan 27, 2022

Two temperature plasma reaction #1166

Merged

Resolve Yaml paper cuts #1179

Merged

ischoegl removed their assignment Jan 28, 2022

This was referenced Jan 30, 2022

Preamble to InterfaceKinetics refactoring #1184

Merged

Rethink handling of 3rd body colliders Cantera/enhancements#133

Closed

speth added 12 commits February 3, 2022 21:57

[Kinetics] Eliminate need for ElementaryReaction3

0959dbf

An elementary reaction is just a Reaction with an ArrheniusRate as its rate coefficient.

[Kinetics] Eliminate need for PlogReaction3

08c44a5

[Kinetics] Eliminate need for BlowersMaselReaction

3e963af

Add a version of warn_deprecated that shows input file location

e60d763

This is useful for cases where the deprecation refers to something that appears in the input file, rather than being a deprecated method.

[Kinetics] Remove '(+M)' from Chebyshev reaction equations

7a0d424

While these reactions are pressure dependent, that dependence is fully incorporated into the rate coefficient and the concentration of third bodies is never directly used in calculating the reaction rate.

[Python] Fix ChebyshevRate.n_temperature and n_pressure

fe76a97

[Kinetics] Eliminate need for ChebyshevReaction3

d10c914

[Python] Rename Reaction.listFromFile to list_from_file

3a47f30

[Python] Rename Species.listFromFile to list_from_file

7511562

[Kinetics] Treat reactions with different rate types as non-duplicates

59617c5

This maintains the standing behavior of treating different reaction types as non-duplicates after making reactions with different rate types use the base Reaction class.

Move Reaction initialization logic out of ReactionFactory

18f6c92

[Python] Make rate a required argument for Reaction constructor

fba2298

[Python] Simplify constructor options for Reaction

9773967

speth force-pushed the simplify-reactions branch 2 times, most recently from 42b6d17 to 2e000fd Compare February 4, 2022 03:58

[Python] Document direct Reaction constructor

8a8a7bb

speth force-pushed the simplify-reactions branch from 2e000fd to 8a8a7bb Compare February 4, 2022 04:40

speth requested a review from ischoegl February 4, 2022 17:53

ischoegl approved these changes Feb 4, 2022

View reviewed changes

ischoegl merged commit 2d772ab into Cantera:main Feb 4, 2022

speth deleted the simplify-reactions branch February 4, 2022 20:50

This was referenced Feb 10, 2022

Yaml2ck script #1185

Closed

ck2yaml crashes when parsing a third body reaction with PLOG #1165

Closed

ischoegl mentioned this pull request Feb 20, 2022

Deprecate redundant kinetics methods in Python API #1202

Merged

5 tasks

ischoegl mentioned this pull request Mar 13, 2022

Eliminate Remaining Reaction Specializations when using MultiRate Evaluators Cantera/enhancements#142

Closed

5 tasks

ischoegl mentioned this pull request May 14, 2022

segfault when using CustomRate #1282

Closed

Uh oh!

Conversation

speth commented Jan 24, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ischoegl commented Jan 25, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

speth commented Jan 25, 2022

Uh oh!

ischoegl commented Jan 25, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov bot commented Jan 26, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ischoegl left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ischoegl Jan 28, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

speth Jan 28, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

ischoegl commented Jan 31, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

speth commented Feb 4, 2022

Uh oh!

ischoegl left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

speth commented Jan 24, 2022 •

edited

Loading

ischoegl commented Jan 25, 2022 •

edited

Loading

ischoegl commented Jan 25, 2022 •

edited

Loading

codecov bot commented Jan 26, 2022 •

edited

Loading

ischoegl left a comment •

edited

Loading

ischoegl Jan 28, 2022 •

edited

Loading

speth Jan 28, 2022 •

edited

Loading

ischoegl commented Jan 31, 2022 •

edited

Loading