Improve output handling in zeroD simulations

[ischoegl]

Fixes #628

Changes proposed in this pull request:

• use CVodesIntegrator to calculate time derivatives via CVodeGetDky
• implement a Taylor series expansion to estimate future states of a ReactorNet object
• specify maximum difference ("atol") for states between advance states
• create new ReactorNet::advanceTowards function as an alternative to ReactorNet::advance, which automatically reduces the time advance when excessive state differences are predicted (time advances are successively reduced by 2^(-1) until differences remain within specified limits. Internal cvode time steps remain unaffected).

All functions are implemented in C++ and are accessible from the Python interface via cython. The new example reactor3.py illustrates differences to the standard ReactorNet example reactor1.py.

[ischoegl]

The fix in commit 38ab0b1 is unrelated to the pull request, and simply fixes a warning during compilation.

[codecov]

Codecov Report

Merging #629 into master will increase coverage by 0.04%.
The diff coverage is 83.33%.

@@            Coverage Diff             @@
##           master     #629      +/-   ##
==========================================
+ Coverage   70.74%   70.78%   +0.04%     
==========================================
  Files         373      373              
  Lines       43339    43447     +108     
==========================================
+ Hits        30660    30754      +94     
- Misses      12679    12693      +14

Impacted Files	Coverage Δ
include/cantera/zeroD/Reactor.h	76.47% <ø> (ø)	⬆️
include/cantera/numerics/CVodesIntegrator.h	0% <ø> (ø)	⬆️
include/cantera/numerics/Integrator.h	3.07% <0%> (-0.32%)	⬇️
src/oneD/Sim1D.cpp	71.72% <100%> (+0.08%)	⬆️
include/cantera/zeroD/ReactorNet.h	78.78% <100%> (+1.36%)	⬆️
src/numerics/CVodesIntegrator.cpp	76.37% <70.58%> (-0.45%)	⬇️
src/zeroD/Reactor.cpp	84.29% <82.6%> (-0.17%)	⬇️
src/zeroD/ReactorNet.cpp	83.68% <94.91%> (+7.34%)	⬆️
src/transport/GasTransport.cpp	90.58% <0%> (+0.2%)	⬆️

---

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 360ac9b...e43bfd5. Read the comment docs.

[ischoegl]

I incorporated first feedback from #628 and renamed some functions for clarity (plus squashed commits). I can add unit tests once I hear back (edit: rudimentary unit tests are now added).

[speth]

I certainly understand the desire for control over the output interval when integrating, and the use of CVODES to extrapolate the solution in order to do so is an intriguing approach. I am a bit concerned about the complexity, however, given that you can get fairly similar control from very a very simple approach such as writing the integration loop as:

tlast = 0
Tlast = 0
while sim.time < 1e-3:
    sim.step()
    if sim.time > tlast + 1e-5 or abs(Tlast - r.T) > delta_T_max:
        tlast = sim.time
        Tlast = r.T
        states.append(r.thermo.state, t=sim.time*1e3)
        print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
                                               r.thermo.P, r.thermo.u))

This doesn't do quite the same thing as the advance_towards method (it results in having output just after either condition is exceeded, rather than just before), but the effect is similar, and I'm not sure that this is any more complicated for the end user than the new set_atol_advance function introduced here.

One particular concern I have is that adding this would make it difficult for us to change to an alternate integrator instead of CVODES at some point in the future, by locking us into requiring an integrator which has the specific feature of being able to access information about the derivatives.

Perhaps what would make this more compelling would be to provide defaults for the tolerances such that decent performance could be obtained without the user having to do anything except replace step with advance_towards.

[bryanwweber]

My skepticism I expressed about this change in the related issue remains. I've been thinking about this, but been swamped with other work lately. I'll try to make some time to collect my thoughts next week. Sorry for the delay! Thanks @speth for poking this back onto my radar.

[ischoegl]

@speth ... thank you for the review, but I will wait for @bryanwweber's comments before continuing.

To briefly restate my motivation for this PR: I specifically wanted to improve the behaviour of advance, and not necessarily do another version of step. The overarching objective was to provide something simple that allows for straight-forward output control. (I do realize that there is the issue of component_name not having a reverse function for ReactorNet objects.)

At the moment, there are three options for output handling:

1. Use step and live with unnecessarily large data files,
2. Use advance at the peril that essential processes may not be properly resolved, or
3. Apply actual python skills to get what you want

Based on experience, (3) is unfortunately not something that can be expected when teaching an UG combustion class at a generic college, (1) is overkill: some students will use Excel to process data (!), leaving (2) with the stated caveats. Hence my desire to improve upon (2).

Regarding the comment about default arguments: it is a reasonable idea. However, would it make sense to merge advance_towards with advance, leaving no step size reduction as the default (i.e. not changing current behaviour unless specified otherwise?). IMHO, introducing a simple way to get to indices within the reactor state vector (i.e. reverse of component_name which remains to be discussed) would create a very intuitive solution.

[ischoegl]

@speth and @bryanwweber ... at this point I resolved all but one of the requests, but am still hoping to obtain feedback. Following up on my comment from a week ago, i.e.

Regarding the comment about default arguments: it is a reasonable idea. However, would it make sense to merge advance_towards with advance, leaving no step size reduction as the default (i.e. not changing current behaviour unless specified otherwise?). IMHO, introducing a simple way to get to indices within the reactor state vector (i.e. reverse of component_name which remains to be discussed) would create a very intuitive solution.

I defined a new method ReactorNet.component_index (currently defined in cython, but can be easily done in C++ also), i.e.

ix = sim.component_index(r.name, 'temperature')

and am still suggesting to merge advance with advance_towards.

[bryanwweber]

Right, thanks for your patience @ischoegl! My biggest concern with this change as implemented is that now we have yet another way for users to advance time in a simulation. Already, I see many people using advance() when they should be using step() and vice-versa. I think the example that @speth showed of the Python necessary to control the output timing (with the caveat that he noted) is easy enough to teach even to a generic US university undergrad. So, I'm not excited about adding a new method to do almost the same thing that we already have two ways of doing.

I'm also not keen on adding additional optional arguments to the advance() function if they will change the current behavior of the function. The current understanding of advance() is that it reaches the end time and changing that behavior by adding some arguments seems like it might cause trouble, i.e., the user has to understand that sometimes this function reaches the end time, and sometimes it doesn't. In this sense, I much prefer the name advance_towards().

That said, the way I would like to see it work if optional arguments are added to advance() is even with tolerances, the integrator still advanced all the way to the passed in time and returned an array with the collected values. Bonus points if it returns a SolutionArray in Python 😄

Overall, I'm a soft 👎 on this change right now. If it can be merged with advance() in a useful way, then I could be more on board. Before this gets merged, a Matlab implementation will also be needed.

On a related but separate note, I think the same task could be accomplished by adding a filter() method to SolutionArray to reduce a set of data given some user-specified tolerances. Could also add a new method or an argument to advance() to automatically return a SolutionArray filled with data taken by stepping.

[bryanwweber]

I want to add some more thoughts to my last comment. I should have started by saying I think this is a valuable idea before launching into my concerns 😄 By "soft" disagreement, I mainly meant that if we could come to consensus that the cognitive overhead for users is offset by the gain in functionality, then I would be on board. My concern is that, as currently implemented, that balance is on the side of too much overhead.

However, I had another thought. I think the general idea of the user saying "I want to integrate for x units" where "x" and "units" are described by the user has great utility. I think framing it in terms of limiting the time that the solver advances is somewhat confusing though, and the idea that the advance_towards() method might not reach the end time confused me for a while when this was first proposed.

So what if the advance() method were modified to take two arguments, a double and a string/int, where the double is the increment and the string/int is a reactor component name/index? The default would remain the same, passing in just the first argument is interpreted as a time, and the solver advances until it reaches that time. But the user could also do advance(50., 'temperature') and advance 50 K in temperature units, and this would have exactly the same conceptual meaning as advance(0.5, 'time') (which is the same as advance(0.5)). That should hopefully avoid some of the concerns I raised in my earlier point while still staying with the spirit of this PR.

This could be largely the same implementation as is currently here, just presented to the user differently. That said, I would be interested in a comparison with an interpolation method rather than (or in addition to) the current gradient estimation method, i.e. similar to the one done by CVode when the passed value of tout is exceeded. Another comparison would be repeatedly calling step() and checking when the limit is exceeded; I'd be interested to know if this would be much slower than CVode with the CV_NORMAL argument if everything was done in the C++ layer (I think that returning control to Python to do the comparison would be quite a bit slower, but I'm curious if leaving everything in C++ would incur much penalty).

[ischoegl]

Thanks for your detailed thoughts, @bryanwweber! I believe the suggestion of advance(50., 'temperature') deserves thought, although time will have to be part of an advance step: otherwise the next output may never be reached. I also agree that including a third way of integrating (beyond step and advance) is not a good idea. C++ allows for different return types as long as the argument list is distinct, which provides an avenue for a solution.

Regarding SolutionArray: I have significant concerns about this approach, as the object is only defined in Cython and lacks a C++ core implementation. Speed of compiled cython code may be adequate (although output handling can be an issue: see e.g. the ic_engine.py example). Nevertheless, there are significant repercussions when it comes to suitable Matlab implementations. I.e. while I am not opposed to the idea of returning SolutionArray objects, there is a lack of infrastructure.

One thing that I am interested in is the acceptance for creating an interface for time derivatives computed by CVODE (which my current implementation hinges on). If there is agreement that it is a useful feature to be broken out in the ReactorNet interface, I'd propose tabling the discussion for advance, and to restrict this PR to implementing the ReactorNet.get_derivative and/or ReactorNet.get_estimate portion.

[bryanwweber]

although time will have to be part of an advance step: otherwise the next output may never be reached

Yes, I avoided some of the difficulty in the implementation of the idea as an initial thought 😃

there are significant repercussions when it comes to suitable Matlab implementations.

The idea would be that C++ would return (a pointer to) an array, and the separate implementations would decide what to do with that, e.g., Cython could turn that into a SolutionArray, Matlab could do whatever was most natural in Matlab, etc.

One thing that I am interested in is the acceptance for creating an interface for time derivatives computed by CVODE

I don't immediately see the utility of a general interface to this beyond the capability that you've shown in this PR, can you provide some other examples?

Once we had converged on an idea, I was actually going to suggest removing at least get_estimate() from the public-facing API, and probably also get_derivative() as well. The names imply that they would be useful for end users to work with, which I don't think is the case. Especially the get_derivative() name would be confusing because there are many kinds of derivatives that would be very useful (e.g., dT/dt, dP/dt, etc.) that this function is not intended to compute (at least, not to my understanding, please correct me if that's wrong).

[ischoegl]

Especially the get_derivative() name would be confusing because there are many kinds of derivatives that would be very useful (e.g., dT/dt, dP/dt, etc.) that this function is not intended to compute

In the current implementation, ReactorNet.get_derivative retrieves the k-th time derivative from CVODE, i.e. it does indeed compute dT/dt for certain reactor types, together with other time derivatives. Specifically, it is the derivative of the composite state vector of the reactor network (which itself is interfaced via ReactorNet.get_state). The function get_estimate would be a Taylor series expansion to estimate the reactor network state at a future time.

PS: one thing I have been wondering about in a loosely related context is the choice of the Feature Scaling approach in advance_to_steady_state. The step size is chosen by CVODE and thus neither constant nor defined? It obviously works, but convergence appears to depend on the CVODE integrator state as much as it does on residual_threshold.

[bryanwweber]

In the current implementation, ReactorNet.get_derivative retrieves the k-th time derivative from CVODE, i.e. it does indeed compute dT/dt for certain reactor types, together with other time derivatives.

OK, so that may have some utility, if a user is interested in dT/dt or d[X]/dt (given the current state vectors).

The function get_estimate would be a Taylor series expansion to estimate the reactor network state at a future time.

Right, so how is that useful for a user? Why wouldn't they just have the integrator do the integration?

advance_to_steady_state()

I'm not sure what you're asking, but the Users' Group is probably a better place for general questions about features (unless you think it should be an Issue on GitHub). I'd like to keep this discussion focused on the features you're proposing. Thanks 😄

[ischoegl]

Alright. I’ll have to think about this for some time: I don’t see a clear path forward and may just limit this PR to get_derivative.

Regarding advance_to_steady_state, it really was a rhetorical question, as convergence is poorly defined in the current code. It’s a minor issue not worth raising by itself, but a more rigorous definition could build on an actual time derivative, which is why I brought it up.

[speth]

[BW] That said, the way I would like to see it work if optional arguments are added to advance() is even with tolerances, the integrator still advanced all the way to the passed in time and returned an array with the collected values. Bonus points if it returns a SolutionArray in Python

I'm not sure there's a useful, general way to implement automatic collection of integrator states. In C++, there's not an obvious, efficient container, and in Python, SolutionArray would only be useful for a network consisting of a single reactor.

[BW] Before this gets merged, a Matlab implementation will also be needed.

Given the current state of the Matlab toolbox, how difficult adding things is, and how often ideas for improving it involve effectively starting over, I would not insist on this.

[BW] That said, I would be interested in a comparison with an interpolation method rather than (or in addition to) the current gradient estimation method, i.e. similar to the one done by CVode when the passed value of tout is exceeded.

That's an interesting thought, and CVODE already provides access to this interpolation via the CVodeGetDky function. Doing something like this would be a little complicated, though. For example, you would have to step past your threshold for (say) temperature change and interpolate back to an output time that satisfies your threshold. However, the next call to advance or advance_towards could end up being either within this interval, or beyond it and require taking additonal steps.

[IS] One thing that I am interested in is the acceptance for creating an interface for time derivatives computed by CVODE (which my current implementation hinges on).

While I can see that the computed derivative function does seem to behave reasonably well in your extrapolation method, I think that this polynomial is only meant to apply to the time interval between the last output time and the current time, given the way the BDF method is defined. So, I'm a little concerned about providing a method whose main purpose is to use this data outside the time interval for which it is applicable.

[bryanwweber]

@speth do you have any thoughts about my concern about adding another user interface method for integration?

However, the next call to advance or advance_towards could end up being either within this interval, or beyond it and require taking additonal steps.

I'm not sure I understand what you're saying here. Why is this a concern?

[ischoegl]

[BW] That said, the way I would like to see it work if optional arguments are added to advance() is even with tolerances, the integrator still advanced all the way to the passed in time and returned an array with the collected values. Bonus points if it returns a SolutionArray in Python

I'm not sure there's a useful, general way to implement automatic collection of integrator states. In C++, there's not an obvious, efficient container, and in Python, SolutionArray would only be useful for a network consisting of a single reactor.

I didn't have a good idea here either, especially as there is no good data structure candidate in C++ implemented at this point (and I am not even sure whether implementing one makes sense: it definitely goes way beyond the original intent of this PR).

[BW] That said, I would be interested in a comparison with an interpolation method rather than (or in addition to) the current gradient estimation method, i.e. similar to the one done by CVode when the passed value of tout is exceeded.

That's an interesting thought, and CVODE already provides access to this interpolation via the CVodeGetDky function. Doing something like this would be a little complicated, though. For example, you would have to step past your threshold for (say) temperature change and interpolate back to an output time that satisfies your threshold. However, the next call to advance or advance_towards could end up being either within this interval, or beyond it and require taking additonal steps.

Agreed. step and advance are calling CVODE with different flags (CV_NORMAL and CV_ONE_STEP, respectively). A solution using interpolation would have to bypass CV_NORMAL completely. One of my questions here is whether this is developing too far away from the philosophy of the underlying integrator.

[IS] One thing that I am interested in is the acceptance for creating an interface for time derivatives computed by CVODE (which my current implementation hinges on).

While I can see that the computed derivative function does seem to behave reasonably well in your extrapolation method, I think that this polynomial is only meant to apply to the time interval between the last output time and the current time, given the way the BDF method is defined. So, I'm a little concerned about providing a method whose main purpose is to use this data outside the time interval for which it is applicable.

Correct. The derivatives are accurate within the time interval of the last CVODE step, although I am not sure how they are internally calculated. What I used for the current implementation was a Taylor series expansion, with all its well known caveats. Time derivatives are, however, also of interest within other contexts: one would be my comment above about the current advance_to_stead_state implementation (which imho should use the derivative instead of an increment over an unspecified time interval in a more rigorous definition).

[ischoegl]

@speth do you have any thoughts about my concern about adding another user interface method for integration?

However, the next call to advance or advance_towards could end up being either within this interval, or beyond it and require taking additonal steps.

I'm not sure I understand what you're saying here. Why is this a concern?

@bryanwweber This is linked to how CVODE is being called (i.e. with flags CV_NORMAL or CV_ONE_STEP). CV_NORMAL uses multiple single steps and interpolates within the last step to the targeted output time. Interpolation is only available for the last step, i.e. the implementation you are suggesting would have to bypass CV_NORMAL entirely. Thus, advance_towards would have to emulate the CV_NORMAL behavior within the C++ layer instead of using the native CVODE code.

[ischoegl]

@speth and @bryanwweber ... I figure that the intent of my original PR was to treat the step size reduction via thresholds (with bounds estimated using a Taylor series expansion). The implementation is straight-forward, and does not bypass CVODE internals.

Summarizing my other thoughts from before, I believe the following Cantera interface makes sense:

• ReactorNet.step() as is
• ReactorNet.advance(tout) as is
• ReactorNet.advance(tout, limit_step=True) replacing advance_towards, with the step thresholds being set by ReactorNet.set_advance_limit (renamed from current ReactorNet.set_atol_advance), or a ReactorNet.advance_limit attribute (renamed from ReactorNet.atol_advance). If no thresholds are set, the behavior defaults to the regular advance behavior.
• ReactorNet.component_index(reactor_name, name) locates components within various ReactorNet structures

That way, not reaching the output time will not be a surprise to the user, which is really the main concern that I can see.

Otherwise, I do not see a simple solution emerging from this discussion.

[speth]

I like this proposal for the interface, using an optional argument to advance.

I think I would prefer the set_advance_limit function to a property, since the most common way to use it would be setting the limit for a specific solution component. This name is definitely an improvement over atol_advance.

While adding the ReactorNet.component_index function would be useful in terms of completing the interface, you could avoid some awkward function calls if the set_advance_limit were also implemented on the Reactor object.

I also agree that having a way of calculating the derivative would be helpful for fixing the flawed implementation of advance_to_steady_state (as part of a separate PR, of course 😄).

[ischoegl]

While adding the ReactorNet.component_index function would be useful in terms of completing the interface, you could avoid some awkward function calls if the set_advance_limit were also implemented on the Reactor object.

This will require some work (and affect more class definitions), but it does make sense.

I also agree that having a way of calculating the derivative would be helpful for fixing the flawed implementation of advance_to_steady_state (as part of a separate PR, of course smile).

I will open an issue (#654)

[ischoegl]

@speth ... the suggested change to the interface is implemented, together with your request to implement set_advance_limit for Reactor objects (as well as for ReactorNet objects). Regarding the comment about properties: I kept the attribute setter methods for the time being (I also kept get_estimate in cython). While some of those methods can be arguably removed from the cython interface, a few are being used internally, i.e. exposing them via cython may make sense.

[ischoegl]

I left a couple of comments.

[ischoegl]

@bryanwweber and @speth ... I believe I took care of everything, but left the two items that were not trivial open (despite considering both as resolved). Let me know whether there is anything else.

[ischoegl]

The last two remaining points are now addressed.

[ischoegl]

@speth ... thanks for the review and comments: I believe everything is addressed. Also let me know if I should rebase the PR (afaik, it's not necessary).