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

Extensible rate: unit handling, serialization, and validation #1478

Merged
merged 26 commits into from Apr 18, 2023

Conversation

speth
Copy link
Member

@speth speth commented Apr 17, 2023

Changes proposed in this pull request

  • Expand handling of AnyMap units in Python
  • Implement ExtensibleRate.get_parameters to enable serialization of user-defined rate types
  • Implement ExtensibleRate.validate to enable validation checks for user-defined types
  • Fix setting of Reaction.rate_units
  • Fix incorrect rate coefficient unit calculations for falloff / chemically-activated reactions
  • Fix confusion where rate coefficient units were being conflated with units of the parameter appearing in the rate expression (specifically, this applies to sticking coefficients)
  • Add tests for rate coefficient units

If applicable, fill in the issue number this pull request is fixing

Contributes to implementing Cantera/enhancements#79

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

From the custom_reactions.py example:

@ct.extension(name="extensible-Arrhenius", data=ExtensibleArrheniusData)
class ExtensibleArrhenius(ct.ExtensibleRate):
    __slots__ = ("A", "b", "Ea_R")
    def set_parameters(self, params, units):
        self.A = params.convert_rate_coeff("A", units)
        self.b = params["b"]
        self.Ea_R = params.convert_activation_energy("Ea", "K")

    def get_parameters(self, params):
        params.set_quantity("A", self.A, self.conversion_units)
        params["b"] = self.b
        params.set_activation_energy("Ea", self.Ea_R, "K")

    def validate(self, equation, soln):
        if self.A < 0:
            raise ValueError(f"Found negative 'A' for reaction {equation}")

    def eval(self, data):
        return self.A * data.T**self.b * exp(-self.Ea_R/data.T)

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
Copy link

codecov bot commented Apr 17, 2023

Codecov Report

Merging #1478 (659df28) into main (6e2a5c0) will increase coverage by 0.08%.
The diff coverage is 90.08%.

❗ Current head 659df28 differs from pull request most recent head 8b27a52. Consider uploading reports for the commit 8b27a52 to get more accurate results

@@            Coverage Diff             @@
##             main    #1478      +/-   ##
==========================================
+ Coverage   70.06%   70.14%   +0.08%     
==========================================
  Files         377      377              
  Lines       57789    58001     +212     
  Branches    19410    19422      +12     
==========================================
+ Hits        40489    40687     +198     
- Misses      14740    14752      +12     
- Partials     2560     2562       +2     
Impacted Files Coverage Δ
include/cantera/base/ExtensionManager.h 27.27% <ø> (ø)
include/cantera/base/Solution.h 87.50% <ø> (ø)
include/cantera/kinetics/ChebyshevRate.h 100.00% <ø> (ø)
interfaces/cython/cantera/preconditioners.pyx 47.36% <ø> (ø)
src/kinetics/PlogRate.cpp 89.52% <ø> (-0.10%) ⬇️
interfaces/cython/cantera/reactor.pyx 90.79% <50.00%> (ø)
interfaces/cython/cantera/speciesthermo.pyx 94.80% <50.00%> (ø)
interfaces/cython/cantera/transport.pyx 72.15% <50.00%> (ø)
src/base/ExtensionManager.cpp 60.86% <50.00%> (-3.84%) ⬇️
include/cantera/kinetics/Arrhenius.h 90.00% <66.66%> (-3.94%) ⬇️
... and 25 more

... and 2 files with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

@speth … changes look good to me; while I did spend a fair amount of time figuring out units when I refactored the rate calculators, your simplifications make sense to me. Good to see the fixes also: for falloff reactions, I was mainly concerned about getting the final units right and making sure that everything matched the CTI implementation; regarding sticking coefficients, they weren’t meant to have units, but I may have adopted/duplicated infrastructure that brought them in the back door.

My main observation in this PR is about the broadening units interface in Python. It is apparent that extensible objects make exposing C++ units a requirement (some bare bones units functionality I had put in place before was merely a narrow interface that allowed for testing). At this point some potential overlaps with pint capabilities start to appear. While the motivations for these unit system interfaces are certainly different, I am wondering about where this would lead eventually? I don’t have a satisfactory answer to this myself at the moment …

PS: from a technical viewpoint, it makes little to no sense to run things through pint for what is proposed here. So my main question is about how the unit handling approaches are differentiated from the user-facing perspective.

@@ -151,8 +153,75 @@ cdef comp_map_to_dict(Composition m):
class CanteraError(RuntimeError):
pass

DimensionalValue = namedtuple('DimensionalValue',
Copy link
Member

Choose a reason for hiding this comment

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

Looks like there is starting to be overlap with what pint can do. I am personally on the fence about where this can/should go, so this is mostly an observation at this point.


cdef class AnyMap(dict):
"""
A key-value store representing objects defined in Cantera's YAML input format.
Copy link
Member

Choose a reason for hiding this comment

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

See above. We potentially have the option to use pint here, instead of grafting it on at a later point?

ischoegl
ischoegl previously approved these changes Apr 18, 2023
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

Apart from the questions I had in my comment, the changes themselves look good to me.

speth added 25 commits April 18, 2023 09:29
Previously, accessing an ExtensibleRate object from a Reaction object
would create a new wrapper object rather than providing access to the
"original" object that is used by the ReactionRateDelegator, and this
wrapper object would not behave the same as the original object except
for methods that passed through via the Delegator.

Now, we return the underlying object, ensuring consistency. This requires
more complex memory management because the ExtensibleRate/Delegator pair
form a reference cycle that can be attached to other objects either from
the C++ or the Python side.
Handling this centrally simplifies the implementation of user-defined
rate functions
The forward rate constants for falloff and chemically activated reactions
already incorporate the third-body dependency, so the units should only
reflect the explicit reactants. This was being compensated for when
setting the units of the high/low pressure rate expressions, so it did
not affect any results.
These differ specifically in the case of sticking reactions, where the
rate coefficient has units like m^3/kmol/s, depending on the reactant
orders while the sticking coefficient itself is dimensionless.
@speth
Copy link
Member Author

speth commented Apr 18, 2023

You're right that there is some significant overlap in what pint can be configured to do and what the C++ UnitSystem class does, but I think there's something to be said for the consistency of always using our internal methods rather than introducing the possibility that we would see differing behavior between Python and C++. I don't anticipate the need for much further expansion of the Python unit handling capability -- this pretty much mirrors what we have in C++, and that has been sufficient for all of our needs on that side.

I think that even if we eventually expand the use of units as in #991 and that becomes a popular way to use Cantera, I'd still expect to retain a lower-level interface that worked in plain mks+kmol values, and you'd probably always want to use that interface for things like ExtensibleRate where having units around would just add a further speed penalty.

I updated the PR to (a) hide the DimensonalValue namedtuple, as this is just an implementation detail (the corresponding C++ Quantity struct is likewise internal) and (b) provide a little more guidance on how the UnitSystem class is meant to be used.

@ischoegl ischoegl merged commit e82f932 into Cantera:main Apr 18, 2023
35 checks passed
@speth speth deleted the extensible-rate-units branch April 18, 2023 14:46
@speth speth mentioned this pull request Apr 24, 2023
4 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants