Eliminate traps associated with accessing ThermoPhase objects in use by Reactors

[speth]

Abstract

For a ThermoPhasein use by a Reactor, the state of the ThermoPhasedoes not necessarily correspond to the state of the Reactor at the current time of the ReactorNet. This can lead to users extracting incorrect state information during integration (see this UG post, for example). In addition, the current requirement that a single ThermoPhasecan be shared across multiple Reactors increases the complexity of evaluating the Reactor's governing equations.

Motivation

After an integration step by ReactorNet, the state of the ThermoPhase object will generally be set to the latest time reached by the CVODES integrator, $t_{int}$. However, this time may be beyond the output time specified by the user (using the advance(t_user) method). If the user accesses the ThermoPhase object directly, they will get the state at $t{int}$, in contrast to getting the state at $t{user}$ by accessing it through the reactor.thermo property, in which CVODES will provide a result interpolated back to the correct time. For simple usage, like plotting the state as a function of time, getting the state at a slightly different time won't noticeably affect results. However, for cases like trying to evaluate sensitivities using a finite difference approach, two reactor networks will not have reached the same value of $t_{int}$, introducing a large amount of noise. This is the phenomenon dubbed "Bjarne's Trap" by @rwest, and has been a source of confusion for a long time.

Because the current structure allows multiple reactors to share a ThermoPhase object, there is a significant amount of complexity in the logic for evaluating the governing equations for all reactors in the network.

The current evaluation logic amounts to the following (in Python pseudocode):

def eval(t, y, ydot):
    # ReactorNet::updateState
    for reactor in all_reactors:
        reactor.thermo.set_TPY(y[...])
        # save the current state information 
        reactor.state = reactor.thermo.state
        # compute auxiliary quantities based on state, for any quantities that may need
        # to be accessed from other reactors, e.g. mass fractions, pressure, or enthalpy
        reactor.local_vars = f(reactor.thermo)

    # ReactorNet::eval
    for reactor in all_reactors:
        reactor.thermo.restore_state(reactor.state)
        ydot[...] = f(reactor.thermo, reactor.kinetics, reactor.local_vars, reactor.neighbors.local_vars)

This complexity has resulted in difficulties for users who want to implement modified governing equations using the ExtensibleReactor framework. For example, see this problem posted on the UG, and this one too.

Possible Solutions

I think the solution to both of these problems to stop having multiple reactors share ThermoPhase/Kinetics objects, and to have ReactorNet.step and ReactorNet.advance automatically synchronize the state of the reactors after each call (by default).

• Thanks to the implementation of AnyMap-based serialization, we can still allow users to use a single Solution object to set up multiple reactors. When initializing the ReactorNet, we can then check to see if any Solution object is being used multiple times, and then clone it to create a distinct Solution for each Reactor.

  • If a gas object is being used only once, don't make a duplicate, and automatically synchronize it after every call to step/advance (TODO: verify if this is necessary for step; not entirely clear when that would be the case)
  • If a gas object is being used for multiple objects, create an internal duplicate for all of the internal uses, leaving the external one corresponding to the initial condition for whichever reactor was set up last. The other reactors will be accessible via the thermo properties of the reactor objects. We will still automatically synchronize these after each step -- not rely on the magic call to syncState when the getter is used.

• The automatic synchronization can be disabled by passing a sync=False option to advance() or step().

• The "simple" alternative is to force this on the user -- we just check and tell the user that they can't use the same Solution object multiple times, but I think the automatic cloning is much more user friendly and shouldn't be that difficult to implement.

After this change, the evaluation logic would look more like:

def eval(t, y, ydot):
    # ReactorNet::updateState
    for reactor in all_reactors:
        # Reactor.updateState
        reactor.thermo.set_TPY(y[...])
        # Calculation of reactor-local variables is optional, for efficiency only

    # ReactorNet::eval
    for reactor in all_reactors:
        # Evaluation can directly access thermo from this reactor and neighboring reactors
        ydot[...] = f(reactor.thermo, reactor.kinetics, reactor.local_vars, reactor.neighbors.thermo)

I'd like to make sure that any changes to how the work of evaluating the reactor governing equations is divided up would provide relatively clean solutions to a few scenarios that have come up on the Users' Group. Namely:

• https://groups.google.com/g/cantera-users/c/0IeoYgEr-Ng: imposing $T(t)$ while keeping other reactor properties (volume, density, etc.) consistent
• https://groups.google.com/g/cantera-users/c/ASrSvruxw-g: using thermo properties in an ExtensibleReactor override override of update_state

[ischoegl]

@speth ... this certainly looks like something we should try to make more user friendly.

Regarding use cases: the first issue (ODE integration) could likely be handled by automatically calling syncState (or something equivalent) at the end? Whenever accessing ReactorBase.thermo, we already have a restoreState in place, so that should be covered?

For the second case, we could also do some thing like

reactor.set_TPY(y[...])

instead of

reactor.thermo.set_TPY(y[...])

and handle things internally. The exact terminology/API is obviously debatable.

At the same time, tackling what amounts to a copy constructor for Solution objects is not necessarily a bad option (would likely be a benefit in other areas also). It's something that is currently missing from the C++ core (although something similar exists for the Python pickle case).

[speth]

Regarding use cases: the first issue (ODE integration) could likely be handled by automatically calling syncState (or something equivalent) at the end?

At the end of what? At the end of step() or advance(), if multiple reactors share the same Solution, which reactor should that Solution's state correspond to?

Whenever accessing ReactorBase.thermo, we already have a restoreState in place, so that should be covered?

This is the main problem: While the thermo property in Python does this synchronization, there's nothing that prevents the user from inadvertently accessing the Solution object directly without using the thermo property.

Don't read too much into my pseudocode structure. I wrote it with some of the functions effectively "inlined" for what I was hoping would improve clarity, but perhaps didn't. I fully expect the manipulation of the ThermoPhase object to be encapsulated within a method of ReactorBase, which will probably still be updateState.

[ischoegl]

This is the main problem: While the thermo property in Python does this synchronization, there's nothing that prevents the user from inadvertently accessing the Solution object directly without using the thermo property.

This is indeed the crux of the problem. From where I stand, it is poor practice if a user has expectations that a shared object should hold a specific state (or attempts to use that implicit pathway to update reactor contents). It is unfortunate that the current design paradigm allows for this. It is getting 'abused' and people report 'unexpected behavior'. Overall, I don't think that this is a path we should encourage.

My concern is that disallowing shared objects in some cases will break existing code (example: it is fairly common to use the same Solution object to initialize several reactors).

As we're already breaking things, we may as well create a clean break and shift to a different paradigm where Solution objects are copied, and cannot be shared (i.e. ReactorBase calls a copy constructor for a new Solution object, rather than attaching the previously created shared object). I would prefer this path as this would not break code where one object is used to initialize multiple reactors (a legitimate use case), while putting an end to bad coding practice (expecting a magic link of Reactor and Solution objects that are not ReactorBase.thermo). This approach would then also apply to 1D objects as well as SolutionArray objects.