Question regarding changing the time step

[radarsat1]

This is not a feature request, just wondering if there is any existing support for changing the time step. I see that there are some old comments in TimeDiscretisation.hpp (from 2004?) that refer to a function setCurrentTimeStep() that doesn't exist any more, however in the current API, I don't see how to do it.

I tried creating a new TimeDiscretisation object and calling Simulation::setTimeDiscretisationPtr(), but that seems to assert on k < _k, so maybe there is a better way to do it?

It's not a critical requirement, it is ok if it is necessary to just re-initialize the entire simulation from scratch. Just that my impression is that this should be possible so I thought it might be worth asking. In any case it seems the documentation in TimeDiscretisation.hpp is out of date so in that sense it is a bug report ;)

[vacary]

The feature was disabled because it was not correctly supported by the simulation classes. Nevertheless, it must be a feature of the simulation classes to be able to change the time step, or to have a variable time--step.

The main difficulty is for linear time invariant systems such as LagrangianLinearTIDS where it assumed that the time--step is constant for computing the iteration matrix W in one step integrator. For the other systems and one--step OSIs, it should be easy to implement a variable time--step, either a) by changing the Timediscretization object, b) using a time-discretization with a non constant time--step. or c) by inserting a new event of type TimeDiscretisationEvent in the stack of events.

I suggest to create examples based on the BouncingBallTS to test these features and to start to reenable this feature.

It will be important for the link of siconos with hybrid modeling languages as Modelica. I can participate to this issue.

To conclude, the better is to avoid to put a boolean that encodes if the time-step is contant of not. It is better to test the value of the current time--step with the previous value.

[radarsat1]

Ok, I agree that having a variable timestep is useful in certain cases, definitely. However for now I was just talking about changing it once, very rarely. For that, it should be sufficient to ignore "correctness" and just change the value in TimeDiscretization, but it does not seem obvious if that's possible. Basically k is then wrong, in the EventManager, if I understand correctly.

[vacary]

Strange again. Normally, we should be able to insert a TimeDiscretizationEvent in the EventManager. We tried to build these classes without referring to a iteration number k. Perhaps, it is not working well.

[vacary]

I just have a look to the EventManager class. It cannot work like that. We have to refactor also this part :-(

[radarsat1]

Ah, shame. In any case, should not be considered high priority. Fix only if it's easy ;) Thanks for taking a look!

[xhub]

I rewrote most of the EventsManager a few years ago, so I should be somehow knowledgeable. The doc is incomplete, but it is a good opportunity to improve that.

I think that @radarsat1's experience with setTimeDiscretisationPtr makes me want to remove that function. It is used nowhere and is not working.

some backgroup on how the EventsManager works:

• there is an event stack. The event at position 0 is the current event and gives the current simulation time.
• When this stack is updated(::update), an update method on the event is used to change the time. If the event time is less than the end of simulation time T, then it is again injected in the stack.

Therefore, it should be possible to implement this functionality by doing the following

• Extend TimeDiscretization to support a variable timestep. This could be done by adding a vector of timesteps, indices and time. The following should old: time[i] = time[i-1] + timestep[i] * index[i], with time[-1] = T0. This scheme allow to reuse the same timestep multiple times. The current time could be computed easily via time = time[end] + timestep * index.
  Another option is to have a callback support that would give the timestep and the next time is then computed using currentTime + oracleTimeStep.

Note that it is super important to minimize the number of arithmetic operation to compute the time. The naive approach time += timestep leads to major drift in both direction because of the IEEE 754. We don't want by default to have this kind of approach. That's why the current way of computing time with a fixed timestep is _t0 + timestep * k.

We could get rid of _k in EventsManager, by having the argument in the update() method being the delta in index instead of an absolute value. We just have to adapt insertEvent(const int type, SP::TimeDiscretisation td), but I can do that. I think this change would be a good enhancement.

[radarsat1]

Thanks for the detailed reply @xhub. I agree that an integer-counter approach is necessary for precision. One solution not requiring a lot of re-engineering could be to simply rescale _k in the event that timestep changes, no?

Yes, setTimeDiscretisationPtr should probably be removed if it cannot be fixed. Indeed it would seem like rather undefined behaviour to swap out the time basis of a Simulation. We could decide to define what happens in that case, or just remove the possibility of undefined behaviour.

I should note that the motivation here is to handle the case in Gazebo where the user uses the GUI to change the physics timestep. In that case it is totally acceptable (to me) to just re-initialize the whole simulation, since it is a rare event. However, for the general case of variable time-step simulations it could be interesting to think about the right way to do this.

(For instance, could a variable-timestep subclass of TimeDiscretization be useful, that uses a floating point accumulator, trading off precision for flexibility?)

[vacary]

Thank you for this detailed discussion.

@xhub, I remember now what was the motivation for this change. Since the drift in the floating point operation is something that can only be measured with a constant time step, I think we can come back to another solution for the general case by doing simply floating point additions.

I agree with @radarsat1 that we have to implement some subclasses of TimeDiscretization to take into account the various cases.

• TimeDiscretization with a constant time step size.
  Let us continue to maintain the rule as it is _t0 + timestep * k to avoid the drift. Should it be the parent class ? It will be easier since we will not have to change all the examples. However, it is curious to define the parent class with the more specific case.
• TimeDiscretization with a highly variable time step.
  In that case the time--step size is changed at nearly all time steps. This is the case with self-adaptive variable time-step size technique in OneStepIntegrator (rule (60) in (http://www.inrialpes.fr/bipop/people/acary/publications/Acary_APNUM2633.pdf)). I do not know in that case what should remain in the TimeDiscretization class since we directly can schedule a event into the stack, but we can design a very
• TimeDiscretization with few resets of the time--step size.
  I think that a simple rescale of _k and a change of the time--step should be sufficient. It can perhaps be implemented directly in an enriched
  TimeDiscretization class with a constant time step size
• TimeDiscretization with a variable known in advance.
  This is the case described by @xhub with the rule time[i] = time[i-1] + timestep[i] * index[i]. I do not know exactly in which cases it can be applied but if it has an interest we can also implement it.

Some remarks on the current implementation of TimeDiscretisation.hpp

• The method

  /** get the timestep \f$t_{k+1} - t_k\f$
    * \param k the index of the timestep
    * \return the current time step
    */
   double currentTimeStep(const unsigned int k);
  

  is somehow misleading. If it is the current time--step that is returned, it should not need k as a parameter.
• The class is only accessed by the EventManager class. So it would be not so difficult to change it. The index _k that registers the index of the current time--step is in the EventManager and not in TimeDiscretization.

I think It would be valuable to create a test in siconos/kernel/src/simulationTools/test/ to test to new features independently of the time integrating a NSDS.

[xhub]

* Vincent Acary <notifications@github.com> [2017-03-07 23:50:12 -0800]:

I agree with @radarsat1 that we have to implement some subclasses of TimeDiscretization to take into account the various cases. - TimeDiscretization with a constant time step size. Let us continue to maintain the rule as it is _t0 + timestep * k to avoid the drift. Should it be the parent class ? It will be easier since we will not have to change all the examples. However, it is curious to define the parent class with the more specific case.

The base class could be TimeDiscretizationBase and be non-instanciable (no public ctor). Then the TimeDiscretization should be named TimeDiscretizationConstant, but that would be ommited by default.

- TimeDiscretization with a highly variable time step. In that case the time--step size is changed at nearly all time steps. This is the case with self-adaptive variable time-step size technique in OneStepIntegrator (rule (60) in (http://www.inrialpes.fr/bipop/people/acary/publications/Acary_APNUM2633.pdf)). I do not know in that case what should remain in the TimeDiscretization class since we directly can schedule a event into the stack, but we can design a very

I think we are missing the end of the sentence :)

- TimeDiscretization with few resets of the time--step size. I think that a simple rescale of _k and a change of the time--step should be sufficient. It can perhaps be implemented directly in an enriched TimeDiscretization class with a constant time step size - TimeDiscretization with a variable known in advance. This is the case described by @xhub with the rule time[i] = time[i-1] + timestep[i] * index[i]. I do not know exactly in which cases it can be applied but if it has an interest we can also implement it.

Yes I don't know if it is of interest either. I was thinking of a case when the timestep does not change too often. There would be a ``standard timstep'' that is used by default. Now I was thinking that if there is a difficult interval to integrate, there could be some external logic to reduce the timestep locally in time so as to go through that situation. But maybe that is not realistic. Thetimestep vector and index vectors are filled during the simulation. I agree that currentTimeStep should not accept any argument. The only place where it is used is in the Control stuff and the code could and should be improved. Feel free to break that.