Codecov Report

Merging #1084 (0766a35) into main (375902b) will increase coverage by 0.00%.
The diff coverage is 66.66%.

@@           Coverage Diff           @@
##             main    #1084   +/-   ##
=======================================
  Coverage   73.45%   73.45%           
=======================================
  Files         362      364    +2     
  Lines       47554    47611   +57     
=======================================
+ Hits        34929    34972   +43     
- Misses      12625    12639   +14     

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 375902b...0766a35. Read the comment docs.

I guess this needs to be expanded to clib/MATLAB 😀 … I’d still appreciate feedback on the implementation approach.

Thanks for the review @decaluwe and thanks for the PR @ischoegl . Personally, I'm not in favor of changing this behavior, for the reasons I mentioned in the issue (even with the legacy flag added). I won't object, though, if @decaluwe or @speth decide to merge

@decaluwe ... thanks for your positive feedback. Regarding literature: Kee just happens to have an equation that is most clear-cut, but other texts are pretty consistent in treating M just like any other species. This implies that it is not commonly understood as being part of the reaction rate constant k. I can point to several examples where this is evident (and am not aware of any references that merge M with k outside of Falloff reactions, which are specifically not the objection here):

• Glassman (3rd ed): Eqs 25/25a
• Kuo: Section 2.3 (Third-Order Reactions)
• Turns (3rd ed): Eqs. 4.24/4.25
• Law: As part of the foundation of Lindemann, Eqs 2.3.3-2.3.6

So, I believe Cantera's inclusion of M in the rate constant of three-body reactions is non-intuitive.

Fwiw, CHEMKIN is consistent with Kee (perhaps not suprisingly): "When a third body is needed, the concentration of the effective third body must appear in the expression for the rate-of-progress variable. ", see Eq. 3-18 (Chemkin 17.0).

Regarding other questions:

Is the idea to switch the default for use_legacy_rate_constants to true, starting with Cantera 2.6?

It would be true for Cantera 2.6 (no change, but warning), and then switched to false after that. For those who want to adopt early and don't want to always switch, there's a scons flag.

Do you intend for the use_legacy_rates functionality to be a permanent feature?

I think the feature can remain, as it's pretty simple to switch. If there's a consensus to remove (after Cantera 2.7/3.0), that would be fine also.

In base/global, you have bool values defined for use_legacy_rate_constants and legacy_rate_constants_used--do I understand correctly that the former is the "set" method, whereas the latter is the "get?"

This is correct. I used the same terminology as the preexisting thermo_warnings feature. It would be simple to change the names.

PS: beyond, I compiled the MATLAB extension in the meantime and fixed a minor glitch.

PPS: just saw that you commented also @bryanwweber ... at this point, I believe it's up to @speth to weigh in. Fwiw, this plays into Cantera/enhancements#87, where the status of M has not been fully resolved.

Thanks, @ingmar - note that I accidentally wrote true in my "switch the default" question (50/50 chance of getting it right 😂 )... but your answer clarified it 100% regardless.

@bryanwweber, can we discuss your reservations a bit more? I'm not really keen on steamrolling over your concerns. To consolidate the convo, you listed:

1. That it was based on 'anecdata' - see @ischoegl additional comments above.
2. Computation speed. As @ischoegl mentioned, getFwdRateConstants merely duplicates a snippet of the updateROP code, rather than offloading calculations from updateROP. I.e. updateROP does not call the change.
3. Code history and breaking any code that depends on the "legacy" way of doing it. Obviously, this is the most serious of the concerns, but as with any compatibility-breaking changes, people will have the deprecation period to make noise about it, if it's going to cause real problems (side note: I always feel like the deprecation message should explicitly mention this as an option). To me, this feels similar to our conversation on standard states, where we are discussing how our convention for a mostly internal routine comports with field conventions. There is "what is the right approach?" and then "If the current approach is not the right one, is it wrong enough to justify the cost of changing it?" To me the costs on this one feel small (but I'll note that I don't use these routines, so my word does not go terribly far, here). I will also note that we are not necessarily "breaking" scripts that have used the current approach. It is also possible that these scripts are using the rate constants as if they are consistent with @ischoegl's definition. In this case, we'd be fixing their scripts 😂

Agree with @ischoegl, though -- I'd like to hear @speth weigh in, since there is still disagreement.

There is "what is the right approach?" and then "If the current approach is not the right one, is it wrong enough to justify the cost of changing it?"

😆 ... to me personally, this always looked like a bug, and - at least to me - it looks wrong enough to be worth fixing; it's not the first time I stumbled across this implementation. Looking over references today confirmed that my recollections square with what I believe the theory says. There may be some other references, but what I listed in my previous post corresponds to my go-to references on my bookshelf. Unfortunately, I don't have Gardiner on hand, but I'd frankly be surprised if it said anything different (lumping the M concentrations together with k makes little intuitive sense to me).

Now, there may be one addition that may facilitate a decision to adopt what I believe is the correct implementation: at the moment, there is no interface to query third-body concentrations outside of C++ (which makes it rather hard to calculate anything based on collision efficiencies). I'd be happy to add this feature to this PR ...

Thanks @decaluwe and @ischoegl. I can provide some more thoughts, as they've evolved over the last few days of thinking about this:

1. Agreed, this does seem to be a standard. Thanks for digging up more references @ischoegl!

2. and 3. I have a few thoughts here still. I personally don't view this as a bug, or even really a problem. The software works as intended and the only downside is that it is somewhat harder to trace the calculations, but at a level that most (for some definition of 'most') people don't care about. As such, I prefer to stick to the "if it ain't broke don't fix it" axiom. I think the appropriate fix would be to document what's happening and why.

   That said, we don't know why the calculation is done this way. Since it goes against convention, I have to assume there was some reason for it when it was first implemented. If that reason is no longer valid, or we're willing to discard that consideration at this point, that's fine. But without knowing why, I'm really hesitant to change it.

One last thought I had in terms of the why is that it may simplify calculations for other phase types where [M] is not a factor. For instance, it may allow you to have a concrete rather than virtual updateROP at the Kinetics level, iff getFwdRateConstants includes [M] in the gas phase. Again, I don't know that this is the case - and that's exactly why I don't feel great about this change 😄 , without knowing why (to the best of our ability to answer that), we can't fully predict what would be unintended consequences. The test suite coverage of some of the other phases is not all that good, so it may miss something like this...

@bryanwweber ... thanks for your additional comments.

I personally don't view this as a bug, or even really a problem. The software works as intended and the only downside is that it is somewhat harder to trace the calculations, but at a level that most (for some definition of 'most') people don't care about. As such, I prefer to stick to the "if it ain't broke don't fix it" axiom. I think the appropriate fix would be to document what's happening and why.

I respectfully disagree here. I believe at this point we can all agree that Cantera output goes against accepted standards. I simply don't see a reason why we should cling to a 'feature' (aka bug) that can be fixed. I believe that the proposed deprecation would take care of it once and for all, instead of assuming that users will read the 'fine print' (which they usually don't). To @decaluwe's point, I believe most people don't even realize that the output is not what they would expect.

Regarding virtual vs concrete: the current implementation of updateROP depends on a few other features that are not visible to Kinetics (falloff, etc.), and is currently defined as a virtual function. The hypothetical case to make it concrete would still rely on other virtual functions, so I believe this is just moving the goal posts. However, this may explain why this historical deviation from the standard was used.

I respectfully disagree here.

Thanks @ischoegl, I'm happy to respectfully agree to disagree, at least until we figure out why this is the way it is (if we ever do) 😄 I certainly see all your points, and I think you're making a strong argument. I'm just personally balancing it a little different way, which is why we have a team to work on this!

Thanks @speth! My concerns are resolved then.

Thank you @speth, @decaluwe and @bryanwweber for your inputs on this, even if we couldn't find complete consensus. Based on @speth's review, I edited all docstrings as well as the warning message, which now reads:

DeprecationWarning: GasKinetics::getFwdRateConstants: Behavior to change after Cantera 2.6;
results will no longer include third-body concentrations for three-body reactions.
To switch to new behavior, use 'cantera.use_legacy_rate_constants(False)' (Python),
'useLegacyRateConstants(0)' (MATLAB), 'Cantera::use_legacy_rate_constants(false)' (C++),
or 'ct_use_legacy_rate_constants(0)' (clib).

Among the interfaces, Fortran is the only one where the switch isn't available: I haven't coded in Fortran since the mid 90's (last millenium!), so I probably shouldn't try.

@bryanwweber ... I hope I hit the mark in my deprecation message and docstring phrasing (I may not be the best word-smith). If you (or anyone else) should have suggestions for improvement, I'd be happy to incorporate.

Thanks @ischoegl the message looks good to me. Will this be produced every time getFwdRateConstants() is called as long as the old behavior is used? If so, it'd be nice to have a way to suppress it but retain the old behavior. One way to do that would be to not print the message if the legacy behavior is explicitly chosen by the user, only showing the message if no option is explicitly specified. Another would be to have a separate option, like suppress_thermo_warnings.

One other comment, and I suppose this applies to all the warn_deprecated() calls in C++, I wonder if it's possible to catch those in Python and turn them into real Python DeprecationWarnings so they can be controlled by Python's warning interfaces? This would be an alternate way to address my comment above, as well, but probably is a lot more work, almost certainly outside the scope of this specific change.

@bryanwweber ... thanks for looking over the messages: this is much appreciated! Per your questions:

Will this be produced every time getFwdRateConstants() is called as long as the old behavior is used?

No, only the first time after each launch of Cantera. For subsequent calls, the deprecation warnings are suppressed.

I suppose this applies to all the warn_deprecated() calls in C++, I wonder if it's possible to catch those in Python and turn them into real Python DeprecationWarnings so they can be controlled by Python's warning interfaces?

That's a really interesting question, and yes, this would go well beyond this PR. It would be great to establish this interface (for any warning), but I don't think that this is a simple endeavor as C++ doesn't come with a warning framework like Python (from what I understand, there are some C++ loggers out there to differentiate exception levels, but that's not the issue I believe). In terms of translating warnings to Python, I am not entirely sure how this would work: as no exceptions are raised that can be caught, any Cython call into C++ would probably have to trigger a query whether there are some buffered logging messages (?) ... I am mostly speculating here as details go well beyond my level of familiarity. @speth has probably more insights here ...

In terms of translating warnings to Python, I am not entirely sure how this would work: as no exceptions are raised that can be caught, any Cython call into C++ would probably have to trigger a query whether there are some buffered logging messages

I think this would require something a bit like how we interface with Python for writing messages to stdout, that is, registering a handler in the Application class and in that handler making whatever the correct calls are to the Python C API. @bryanwweber, you're always welcome to write up an enhancement proposal 😁.

@speth - thanks for the feedback/review!

@bryanwweber, you're always welcome to write up an enhancement proposal 😁.

I actually pre-empted this by a discussion here: Cantera/enhancements#112. Good to hear that there may be an easier route!