Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reaction rate factories #1061

Merged
merged 19 commits into from
Aug 17, 2021
Merged

Reaction rate factories #1061

merged 19 commits into from
Aug 17, 2021

Conversation

ischoegl
Copy link
Member

@ischoegl ischoegl commented Jun 14, 2021

Changes proposed in this pull request

Create factory constructors in anticipation of refactoring of Falloff reactions and corresponding rates

  • Implement C++ ReactionRateBase::newMultiRate instantiation
  • Implement C++ newRate factory constructor
  • Generalize Python ReactionRate interface

The Interface now mostly replicates what already exists for Reaction objects. Addresses Cantera/enhancements#87.

Closes #1071 (associated discussion)

Example

YAML input for ReactionRate uses the same information as Reaction. Non-rate specific information is ignored, and rate types are recognized using aliases.

In [1]: import cantera as ct
   ...:  gas = ct.Solution('kineticsfromscratch.yaml')

In [2]: yaml = """
   ...:     equation: H2 + O2 <=> 2 OH
   ...:     type: pressure-dependent-Arrhenius
   ...:     rate-constants:
   ...:     - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04 cal/mol}
   ...:     - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04 cal/mol}
   ...:     - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04 cal/mol}
   ...:     - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04 cal/mol}
   ...:     """

In [3]: rr = ct.ReactionRate.fromYaml(yaml, gas)
   ...: rr
Out[3]: <PlogRate at 7f73992fa060>

In [4]: dd = rr.input_data
   ...: dd
Out[4]: 
{'rate-constants': [{'A': 1.2124e+16,
   'Ea': 45491376.80000001,
   'P': 1013.2500000000003,
   'b': -0.5779},
  {'A': 4.9108e+31, 'Ea': 103649395.2, 'P': 101324.99999999999, 'b': -4.8507},
  {'A': 1.2866e+47, 'Ea': 166508556.0, 'P': 1013249.9999999992, 'b': -9.0246},
  {'A': 5.9632e+56, 'Ea': 220076726.4, 'P': 10132499.999999985, 'b': -11.529}],
 'type': 'PlogRate'}

In [5]: ct.ReactionRate.from_dict(dd, gas)
Out[5]: <PlogRate at 7f73ce3d7660>

For the time being, I am sticking with the usual nomenclature (see #1053)

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 ischoegl marked this pull request as draft June 14, 2021 21:18
@codecov
Copy link

codecov bot commented Jun 14, 2021

Codecov Report

Merging #1061 (39d4577) into main (3dc940d) will increase coverage by 0.51%.
The diff coverage is 92.55%.

Impacted file tree graph

@@            Coverage Diff             @@
##             main    #1061      +/-   ##
==========================================
+ Coverage   72.92%   73.44%   +0.51%     
==========================================
  Files         357      364       +7     
  Lines       47120    47595     +475     
==========================================
+ Hits        34363    34955     +592     
+ Misses      12757    12640     -117     
Impacted Files Coverage Δ
include/cantera/kinetics/MultiRate.h 82.60% <ø> (ø)
include/cantera/kinetics/Reaction.h 100.00% <ø> (ø)
src/kinetics/Reaction.cpp 89.20% <85.71%> (-1.24%) ⬇️
src/kinetics/ReactionRateFactory.cpp 93.33% <93.33%> (ø)
include/cantera/kinetics/ReactionRate.h 61.33% <100.00%> (-0.51%) ⬇️
include/cantera/kinetics/ReactionRateFactory.h 100.00% <100.00%> (ø)
src/kinetics/BulkKinetics.cpp 91.00% <100.00%> (-0.08%) ⬇️
src/kinetics/ReactionRate.cpp 94.25% <100.00%> (-1.75%) ⬇️
src/kinetics/RxnRates.cpp 92.66% <100.00%> (+0.06%) ⬆️
src/base/global.cpp 72.38% <0.00%> (-2.86%) ⬇️
... and 44 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 3dc940d...39d4577. Read the comment docs.

@speth
Copy link
Member

speth commented Jun 15, 2021

We can avoid the need for the MultiRateFactory class by adding a member function to the ReactionRate class, as suggested in my original outline for the MultiRate structure in #982. Namely, for each type of ReactionRate, all that's needed is a function:

virtual shared_ptr<MultiRateBase> newMultiRate() const {
    return make_shared<MultiBulkRate<ArrheniusRate, ArrheniusData>>();
}

Then, the replacement else-if tree in addReaction would work by just calling rate->newMultiRate().

@ischoegl ischoegl marked this pull request as ready for review June 15, 2021 21:00
@ischoegl
Copy link
Member Author

@speth ... thanks for the suggestion - this is indeed simpler.

Copy link
Member

@speth speth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, @ischoegl, most of this looks good to me. I think the only substantial concern I have is with figuring out the rate constant units in the newRate(AnyMap, Kinetics) function.

src/kinetics/RateFactory.cpp Outdated Show resolved Hide resolved
include/cantera/kinetics/RateFactory.h Outdated Show resolved Hide resolved
src/kinetics/RateFactory.cpp Outdated Show resolved Hide resolved
src/kinetics/ReactionRate.cpp Outdated Show resolved Hide resolved
src/kinetics/Reaction.cpp Outdated Show resolved Hide resolved
interfaces/cython/cantera/reaction.pyx Outdated Show resolved Hide resolved
interfaces/cython/cantera/reaction.pyx Outdated Show resolved Hide resolved
interfaces/cython/cantera/reaction.pyx Outdated Show resolved Hide resolved
interfaces/cython/cantera/reaction.pyx Outdated Show resolved Hide resolved
src/kinetics/RateFactory.cpp Outdated Show resolved Hide resolved
Copy link
Member

@bryanwweber bryanwweber left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few small style comments, thanks @ischoegl

interfaces/cython/cantera/reaction.pyx Outdated Show resolved Hide resolved
interfaces/cython/cantera/test/test_reaction.py Outdated Show resolved Hide resolved
@ischoegl
Copy link
Member Author

ischoegl commented Jul 15, 2021

@speth and @bryanwweber - thank you for your reviews ... I addressed / responded to all comments at this point.

Regarding the unit controversy, I believe that there are some inconsistencies that go beyond this PR (see new issue #1071) - the constructor in question is useful for the creation of Reaction objects from scratch, where unit handling is not clear.

@ischoegl ischoegl mentioned this pull request Jul 15, 2021
11 tasks
@bryanwweber bryanwweber dismissed their stale review July 16, 2021 01:14

Thanks @ischoegl, my comments have been address so I'm dismissing my review. @speth can approve since his comments were more about the implementation.

@ischoegl
Copy link
Member Author

Thanks, @bryanwweber!

@ischoegl
Copy link
Member Author

ischoegl commented Jul 17, 2021

@speth ... I believe I resolved the issue with units by removing the Kinetics parameter. Regarding the question about where the creation of a stand-alone ReactionRate object could be useful, I believe it is critical for the creation of a Reaction object from scratch, e.g.

In [1]: import cantera as ct
   ...: 
   ...: reactants = {'O': 1., 'H2': 1}
   ...: products = {'H': 1., 'OH': 1.}
   ...: rxn = ct.ElementaryReaction(reactants, products)
   ...: rxn
   ...: 
Out[1]: <ElementaryReaction: H2 + O <=> H + OH>

In [2]: rxn.rate # reaction rate is not set up (yet)
Out[2]: <ArrheniusRate at 7f38eb99ef00>

In [3]: rxn.rate(300.)
Out[3]: nan

In [4]: rxn.rate = ct.ReactionRate.from_yaml("rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol}")
   ...: rxn.rate
   ...: 
Out[4]: <ArrheniusRate at 7f38eb99ef60>

In [5]: rxn.rate(300.)
Out[5]: 5195.44737317942

@ischoegl ischoegl removed the request for review from bryanwweber July 17, 2021 17:05
@ischoegl ischoegl force-pushed the reaction-rate-factories branch 2 times, most recently from b689e27 to f4595d2 Compare July 22, 2021 14:38
@ischoegl
Copy link
Member Author

ischoegl commented Jul 22, 2021

@speth ... the last commit blocks an edge case where a user would specify a unit system that would produce rate_units that have a factor other than 1.0 (which is probably the main take-away of the discussion in #1071):

In [2]: rxn.rate = ct.ReactionRate.from_yaml("""
   ...: rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol}
   ...: units: {length: cm, quantity: mol}
   ...: """)

[...]

CanteraError: 
*******************************************************************************
InputFileError thrown by ReactionRateFactory::newReactionRate:
Error on line 2 of input string:
Alternative units for 'length' or 'quantity` are not supported when creating
a standalone 'ReactionRate' object.
|  Line |
|     1 | 
|     2 | rate-constant: {A: 38.7 kmol/cm^3, b: 2.7, Ea: 6260.0 cal/mol}
>     3 > units: {length: cm, quantity: mol}
                 ^
*******************************************************************************

With this update, I believe that all comments are taken care of.

PS: there is an exception though:

ct.ReactionRate.from_yaml("rate-constant: {A: 4.0e+21 cm^6/mol^2/s, b: 0.0, Ea: 1207.72688}")

which is a bit harder to catch (see units-custom.yaml).

@ischoegl ischoegl force-pushed the reaction-rate-factories branch 3 times, most recently from 328be0d to f41f20d Compare July 22, 2021 17:17
For stand-alone ReactionRate objects, the pre-exponential factor needs to be
consistent with default Cantera units (length in 'm' and quantity in 'kmol').
In other words, it must not be specified by explicit units or an incompatible
unit system.
@ischoegl
Copy link
Member Author

ischoegl commented Jul 22, 2021

Not too hard after all. I believe this is done now ...

In [2]: rate = ct.ReactionRate.from_yaml("rate-constant: {A: 4.0e+21 cm^6/mol^2/s, b: 0.0, Ea: 1207.72688}")
---------------------------------------------------------------------------
CanteraError                              Traceback (most recent call last)
<ipython-input-2-d77e48358ab6> in <module>()
----> 1 rate = ct.ReactionRate.from_yaml("rate-constant: {A: 4.0e+21 cm^6/mol^2/s, b: 0.0, Ea: 1207.72688}")

/src/interfaces/cython/cantera/reaction.pyx in cantera._cantera.ReactionRate.from_yaml()
    109         cdef CxxAnyMap any_map
    110         any_map = AnyMapFromYamlString(stringify(text))
--> 111         cxx_rate = CxxNewReactionRate(any_map)
    112         return ReactionRate.wrap(cxx_rate)
    113 

CanteraError: 
*******************************************************************************
InputFileError thrown by Arrhenius::setParameters:
Error on line 1 of input string:
Specification of units is not supported for pre-exponential factor when
creating a standalone 'ReactionRate' object.
|  Line |
>     1 > rate-constant: {A: 4.0e+21 cm^6/mol^2/s, b: 0.0, Ea: 1207.72688}
                                   ^
*******************************************************************************

Copy link
Member

@speth speth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the updates to this, @ischoegl. After the extended discussion on #1071, I think most of the complicated issues here have already been worked out without introducing too much new complexity. My remaining comments are mostly places on things that should either simplify the current implementation slightly or make some further generalizations that I suspect we both have in mind easier.

src/kinetics/Reaction.cpp Outdated Show resolved Hide resolved
src/kinetics/Reaction.cpp Show resolved Hide resolved
src/kinetics/Reaction.cpp Outdated Show resolved Hide resolved
src/kinetics/ReactionRateFactory.cpp Outdated Show resolved Hide resolved
src/kinetics/ReactionRateFactory.cpp Outdated Show resolved Hide resolved
src/kinetics/ReactionRateFactory.cpp Outdated Show resolved Hide resolved
@ischoegl
Copy link
Member Author

@speth ... thank you for your re-review and the interesting comments. I ended up adopting all your suggested changes.

[...] one of the key benefits of the generic ReactionRate class is breaking this requirement that there is a specific rate parameterization that goes along with a given Reaction class.

I honestly was still thinking along the old lines too much here. I certainly realized that the revised Reaction3 specializations were pretty much hollowed out (i.e. not many differences), but following that thought it will become possible to at least partially ditch having specialized Reaction classes altogether. As an example, there's nothing wrong with an elementary reaction having a Chebyshev rate object. This opens up a lot of interesting possibilities and simplifications ...

@ischoegl ischoegl requested review from speth and removed request for speth August 12, 2021 15:44
Setting the units factor to 0 is a simpler solution to detect ReactionRate
objects that do not have an associated Reaction defined (yet).

Not calling the default constructor for ElementaryReaction3 from the
constructor avoids the creation of rate objects that are immediately
overwritten.
Removing checks for consistency between (current) Reaction3 objects
and ReactionRate objects makes the framework more flexible.
As ReactionRate objects have consistent type names, specializations for
ElementaryReaction3::getParameters, PlogReaction3::getParameters and
ChebyshevReaction3::getParameters are no longer needed.
@ischoegl
Copy link
Member Author

ischoegl commented Aug 12, 2021

Realized that shifting the emphasis away from Reaction (using consistent naming) makes specialized getParameters mostly unnecessary (except for ThreeBodyReaction3, at least for the moment). Which means that ElementaryReaction3, PlogReaction3 and ChebyshevReaction3 are essentially the same (except for minor differences in constructors). Not bad ... 🎉

I believe this is ready to merge now. Getting rid of some Reaction specializations (while possible) is likely a longer discussion for another PR.

PS: Changing the type string output created a minor paper cut, which is also fixed.

@ischoegl ischoegl requested a review from speth August 12, 2021 16:44
Copy link
Member

@speth speth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, @ischoegl, this looks good to me. I had just one very minor nitpick, commented on below.

Regarding the fate of the Reaction class, I expect that the only distinction ultimately worth keeping will be between bulk phase reactions and interface/electrochemical reactions. But even then, some of the extra parameters and functionality may fit better within the ReactionRate classes.

src/kinetics/ReactionRate.cpp Outdated Show resolved Hide resolved
@ischoegl
Copy link
Member Author

@speth ... thanks for the reviews! I removed the redundant returns with the last commit.

Regarding the fate of the Reaction class, I expect that the only distinction ultimately worth keeping will be between bulk phase reactions and interface/electrochemical reactions. But even then, some of the extra parameters and functionality may fit better within the ReactionRate classes.

👍 ... I see it the same way - the YAML type field will essentially be pushed through to the ReactionRate object, which should be able to take care of (almost) everything.

@speth speth merged commit ec67ed1 into Cantera:main Aug 17, 2021
@ischoegl ischoegl deleted the reaction-rate-factories branch August 17, 2021 20:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Inconsistency in ReactionRateCoeff units calculation
3 participants