Solver/integrator substepping

[RossNordby]

As of now, the v2 timestep loop looks like:

for (0 ... infinity) 
{ 
   position integration
   velocity integration
   bounding box update
   collision detection
   solve constraints
}

To increase stability, it's often recommended to use more timesteps of a shorter duration if you can afford it. Every stage in the above is executed in each timestep. However, it would be relatively easy to shuffle things around:

for (0 ... infinity) 
{ 
   bounding box update
   collision detection
   for (0 ... N) 
   {  
      solve constraints
      position integration
      velocity integration
   }
}

where each substep covers a duration equal to dt / N. The reason for doing this is to focus the extra work of additional time steps where it is most beneficial- constraints. In a simulation where collision detection and constraint solving take roughly the same amount of time, taking 4 substeps for constraints will only be about 2.5x more expensive as a whole compared to 4x more expensive for taking 4 full timesteps. In other words, this is just amortizing the cost of collision detection (and the other minor stages) over multiple integration steps.

Also, increasing timestep frequency provides a much larger improvement in quality per CPU cycle than increasing solver iterations, so a simulation using solver substepping would likely be able to use fewer velocity iterations. In many cases, substepping would increase cost sublinearly, and (compared to using a huge number of iterations) it can even be cheaper.

There are many possible permutations each with their own tradeoffs. The tradeoffs belong to one of two categories:

1. The choice of order exposes a different part of the loop to a user working primarily in between timesteps. For example, if position updates occur first (as they do in the current timestep layout), then changing velocities outside of the timestep results in a direct integration of those velocities without interference from the solver. This can force the user to instead work within callbacks.
2. Performance. For example. splitting position and velocity integration from each other or from bounding box updates results in higher overall memory bandwidth requirements because the body properties get iterated over multiple times and L1/L2 cache will get evicted in between runs.

Remaining questions:

1. What is the ideal permutation of stage execution for usability? Is having contacts that match the post-timestep visible positions more valuable than being able to modify velocity outside of the timestep without worrying about it being directly integrated into position?
2. How do you minimize the cost of bounding box updates? They only need to occur once per narrow phase invocation, but they depend upon a velocity estimate spanning the full timestep.
3. Should contact data be updated after each substep? Without updates to match position integration, the penetration depth will remain the same even if the objects have moved toward or away from each other along the contact normal. This would not be catastrophic (or even meaningfully worse than not substepping to begin with), but the quality improvement may be enough to warrant the cost of incrementally updating contact depths.
4. Is there any value in a variable number of substeps? For example, using a time accumulator to trigger extra solver substeps separately from narrow phase updates. Leaning towards no due to the significant computation variability and complexity.
5. Is it worth making execution order configurable? Completely freeform order isn't practical due to some complicated dependencies, but it wouldn't be that much work to have four or five standard modes of operation. The most difficult part would be exposing this to the user in a way that isn't super gross.

[RossNordby]

Added in b1f3725c3d802c250fb33c7eadc4f51481ec1883. Question of ideal stage ordering is left to the user; ITimestepper interface abstracts the choice.