Update three body logic

[ischoegl]

Changes proposed in this pull request

• Fix cases where three-body detection logic introduced in Remove reaction specializations #1333 fails

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

In [2]: ct.Reaction(
   ...:     equation="CH2OCH + M <=> CH2CHO + M", rate={"A": 5.0e09, "b": 0.0, "Ea": 0.0}
   ...: )
Out[2]: CH2OCH + M <=> CH2CHO + M    <three-body-Arrhenius>

plus several edge cases where three-body reactions are explicitly or implicitly defined, e.g.

In [3]: gas = ct.Solution("h2o2.yaml")
   ...: _yaml = """
   ...:         equation: H2 + O2 <=> H2 + O2
   ...:         type: three-body-Blowers-Masel # or type: Blowers-Masel
   ...:         rate-constant: {A: 38700, b: 2.7, Ea0: 2.619184e4 cal/mol, w: 4.184e9 cal/mol}
   ...:         efficiencies: {O2: 1}
   ...:         default-efficiency: 0.
   ...:         """
   ...: ct.Reaction.from_yaml(_yaml, gas)
Out[3]: H2 + O2 <=> H2 + O2    <three-body-Blowers-Masel>

or

In [4]: rxn = ct.Reaction(
   ...:     equation="H2 + O2 <=> H2 + O2",
   ...:     rate={"A": 5.0e09, "b": 0.0, "Ea": 0.0},
   ...:     third_body=ct.ThirdBody("O2"),
   ...: )
   ...: rxn
Out[4]: H2 + O2 <=> H2 + O2    <three-body-Arrhenius>

In [5]: rxn.input_data
Out[5]:
{'equation': 'H2 + O2 <=> H2 + O2',
 'rate-constant': {'A': 5000000000.0, 'b': 0.0, 'Ea': 0.0},
 'efficiencies': {'O2': 1.0},
 'default-efficiency': 0.0}

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

[codecov]

Codecov Report

Merging #1338 (ced29d7) into main (d60e59e) will increase coverage by 0.10%.
The diff coverage is 83.58%.

@@            Coverage Diff             @@
##             main    #1338      +/-   ##
==========================================
+ Coverage   67.93%   68.03%   +0.10%     
==========================================
  Files         316      327      +11     
  Lines       41948    42652     +704     
  Branches    16853    17180     +327     
==========================================
+ Hits        28498    29020     +522     
- Misses      11186    11345     +159     
- Partials     2264     2287      +23     

Impacted Files	Coverage Δ
include/cantera/kinetics/Falloff.h	81.14% <ø> (ø)
include/cantera/kinetics/Reaction.h	80.00% <ø> (ø)
src/kinetics/Reaction.cpp	80.62% <83.58%> (+0.53%)	⬆️
include/cantera/numerics/FuncEval.h	57.89% <0.00%> (-26.73%)	⬇️
src/numerics/FuncEval.cpp	29.16% <0.00%> (-20.84%)	⬇️
include/cantera/zeroD/Reactor.h	59.09% <0.00%> (-5.91%)	⬇️
include/cantera/zeroD/ReactorBase.h	58.53% <0.00%> (-0.44%)	⬇️
src/zeroD/ReactorFactory.cpp	100.00% <0.00%> (ø)
include/cantera/thermo/Phase.h	83.67% <0.00%> (ø)
include/cantera/zeroD/ReactorNet.h	76.66% <0.00%> (ø)
... and 19 more

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

[ischoegl]

I believe this is ready for a review.

[bryanwweber]

I have one very minor suggestion, which may not be an improvement. I also wonder if you can clarify, in the docstring and the error message, what "ambiguous" means in this context. What should the user do to correct this problem?

[ischoegl]

@bryanwweber ... thanks! I updated the comments and error messages. Let me know if this clarifies things sufficiently.

[ischoegl]

PS: the only question I have is about another edge case. At the moment, the logic requires either three reactants or three products to detect a three-body reaction and thus misses the following

In [6]: ct.Reaction(
   ...:     equation="CH2OCH + O2 <=> CH2CHO + O2", rate={"A": 5.0e09, "b": 0.0, "Ea": 0.0}
   ...: )
Out[6]: CH2OCH + O2 <=> CH2CHO + O2    <Arrhenius>

... it is certainly perceivable that this should be a three-body-Arrhenius instead? Fwiw, the case should be equivalent to

In [7]: rxn = ct.Reaction(
   ...:     equation="CH2OCH + M <=> CH2CHO + M",
   ...:     rate={"A": 5.0e09, "b": 0.0, "Ea": 0.0},
   ...:     third_body=ct.ThirdBody("O2"),
   ...: )
   ...: rxn
Out[7]: CH2OCH + M <=> CH2CHO + M    <three-body-Arrhenius>

In [8]: rxn.input_data
Out[8]:
{'equation': 'CH2OCH + M <=> CH2CHO + M',
 'rate-constant': {'A': 5000000000.0, 'b': 0.0, 'Ea': 0.0},
 'efficiencies': {'O2': 1.0},
 'default-efficiency': 0.0}

PPS: regardless, serialization works and I'd rather not introduce another set of cases where duplicates are created where there weren't any before.

[ischoegl]

I think this is finally ready. At this point I believe that unit tests cover everything I can think of - once this is merged, we should be good to bump the alpha version with #1337

[speth]

Thanks for working on these refinements, @ischoegl.

Beside the suggestions below, I was thinking that we should probably limit the use of the efficiencies and default-efficiency fields to the case where third body is M, to avoid allowing the rather confusing example below:

equation: H + O2 + N2 = HO2 + N2
rate-constant: {A: 1, b: 1, Ea: 200}
efficiencies: {N2: 10.0}
default-efficiency: 1

which is currently allowed and seems to be equivalent to writing the same reaction with M as the third body.

[ischoegl]

@speth ... thanks for the review! I adopted recommendations and tweaked the logic so that your example throws the following error:

*******************************************************************************
InputFileError thrown by ThirdBody::setParameters:
Error on line 5 of input string:
Invalid default efficiency for explicit collider N2;
value is optional and/or needs to be zero
|  Line |
|     1 |
|     2 | equation: H + O2 + N2 = HO2 + N2
|     3 | rate-constant: {A: 1, b: 1, Ea: 200}
|     4 | efficiencies: {N2: 10.0}
>     5 > default-efficiency: 1
                              ^
*******************************************************************************

(it does make sense that the default efficiency has to be zero in this case).

On the other hand, I'd like to keep the ability to change the efficiency for N2. Among others, specifying the efficiency does allow for the disambiguation of the following:

equation: CH2 + O2 <=> CH2 + O2
rate-constant: {A: 5.0e+9, b: 0.0, Ea: 0.0}
efficiencies: {O2: 1.}

I.e. this is handled as a three-body reaction, where the user specifies which of the two participating species is treated as a third body. (As an aside: I realize that the reaction makes little sense; I had originally blocked those but the feedback was that these reactions should be allowed.)

[speth]

Thanks, @ischoegl. This looks good to me.