I would be interested in having a doStep() method implemented by ODEIntegrator and FieldODEIntegrator to have more flexibility in the control of the integration flow.
This method should take as input a (Field)ODEState or (Field)ODEStateAndDerivative and a time step h and return the (Field)ODEStateAndDerivative after one integration step.
For variable step size integrators, the output should be the state at the end of the first accepted step, and the method should also return a guess for the next time step h.
As reference, this concept is implemented for example in boost::numeric::odeint, see Steppers - 1.80.0
As an alternative, having a (Field)ODEStepHandler capable to reset the current state (as (Field)ODEEventHandlers can do) should fit my use case too, but might be less generic than the above solution.
This is not easy.
The features that make this difficult are:
adaptive stepsize integrators
multistep integrators
events detection
Adaptive stepsize integrators compute the step by themselves, so what should they do when the step specified by caller is different from the one they would use? Should they consider the user step as mandatory and use it? Should they consider the user step as a target time and perform if needed several substeps with their own step size before reaching this target?
Multistep integrators need several prior steps to advance, so how should they work in single-stepping calls? Should they consider it is a restart and build up the history internally (this is how they start right now)? Should they also preserve history from one call to the next? What about the case user asks for different steps with varying durations in row?
Events detection works by preserving some internal state between steps (the previous values of the event defining switching functions). How should this internal state data be preserved?
As you see, all the problems above are related to a stateful representation used within the integrator. This information is stored and managed properly for now because the integrate method is called once so the integrator knows how to initialize first, then loop and manage the state between each iteration in the loop. IF we break this loop, we have to pass this information somehow and be sure than one iteration will be consistent with the previous one even if user decided to change everything in between (including introducing discontinuities, or recomputing the same step, or event going back in time).
My use case requires to abandon the while loop over time steps implemented in ODEIntegrator#integrate() and replace it by a custom one, defined on a per-application basis, which would look like:
public ODEStateAndDerivative integrate(ODEState y0, double tf) {
// some initialization ...
double h = ... // compute initial time step
ODEStateAndDerivative y = y0;
do {
ODEStateAndDerivative y = doStep(y, h);
// h should be updated to contain the next time step guessed by variable step size integrators
y = ...; // perform some manipulation on y
if (...) { // check some condition on y
break;
}
} while (y.getTime() < tf);
return y;
}
Concerning adaptive step size integrators, doStep(y,h) should do the following:
compute y1 from y,h
estimate the error
a. if err<tol compute h_new for the next step and return y1,h_new
b. else reject y1 and reduce h, then come back to 1.
I.e., I am interested in extracting the state at the end of the first successful step, as well as a guess for the next time step.
That’s why I think the same goal could be achieved by implementing a resetState method in ODEStepHandler similar to the one in ODEEventHandler.
I am not familiar with multi-step integrators, so I won’t comment on that.
I think event detection is meaningless in this case.
Finally, I am more interested in the Field version of all this rather than the double one.
Just a random thought: couldn’t you just call integrate with a target time at current time + h?
If you use Runge-Kutta based integrator (like DormandPrince853), the overhead at integration start is small, an even negligible if you don’t have many event detectors and step handlers.
For multi-step methods (like Adams-Bashforth or Adams-Moulton), the starting overhead is much larger.
For fix step size methods yes, and I am already doing that.
But, since I need to access the state after each step (independently from the size of h), variable step size integrators would lose their advantage. I would in fact be forced to set h to the largest value that guarantees the accuracy in each step (to ensure just one step is performed at each integrate method call), which can be considerably smaller that the adaptive step size in regions where the dynamics are slower. I don’t think that having an external algorithm to adapt the next target time to restore some form of adaptativity is a good idea.
For now I am using only RK-type methods, not multi-step ones.
So in the adaptive stepsize case, wouldn’t it be more logical to have a doStep(y) method without step size?
I am leaning towards a solution based on step handlers rather than having a doStep method, it seems more generic to me but I still have to figure out how to do that.
Another thing we have to consider is how it interacts with events. As an example, what happens if you change the state y and it changes the sign of an event detector ? Should we generate an event at the beginning of next step or ignore the sign change and consider we have done something similar to Action.RESET_STATE?
If the integrator stores h internally then yes, doStep(y) with just y as input is more logical.
I agree with a solution based on step handlers, in which case I would do the following:
handleStep() should return an Action based on the input ODEStateInterpolator, just like event detectors do
a. Action.CONTINUE mimics the current behaviour, i.e. we can log the state or some quantities derived from it but not influence the integration flow
b. Action.STOP stops the integration
c. Action.RESET_STATE resets the ODEState before resuming the integration
d. possibly others?
ODEStepHandler should implement a method resetState which is called when (c) is returned
With (a) and (b) there is no conflict with event detectors (the first is “passive”, the second just stops the process).
For (c) I would act as if an event that throws Action.RESET_STATE occurred exactly at the end of the step, so ignore the sign change.