Particle System Code - Is Verlet and Midpoint implemented correctly?

I’ve just begun familiarising myself with the blender codebase, in particular the particle system, and I noticed something odd.

Intro for the unfamiliar - skip if you know how numerical integrators work:

So with separable symplectic (Hamiltonian/Energy conserving) integrators, you have two integrator operators: Kick and Drift. Normally you combine these in a weighted fashion to integrate the system. ie. Kick-Drift-Kick with a timestep weight of 0.5-1.0-0.5 gives you the Leapfrog Method (side note: blender uses Leapfrog which has been incorrectly named Midpoint).

Verlet uses a Drift-Kick method, but updates the velocity in the Kick step using the average acceleration from before and after the Drift step.

But Blender doesn’t do this, or at least from my understanding of the code, it does not do this at all.

The issue here being that the particle system treats each particle as separate from the other particles and applies the operators to each particle. This isn’t correct as the acceleration after the Drift operation depends on not only the current particles position, but also all the other particles.

Let’s take the leapfrog method as an illustrative example with its Kick-Drift-Kick pattern.

The current strategy is to apply the Kick-Drift-Kick operators to a particle and calculate an updated acceleration each time one of these steps is applied. Then to repeat this for all the particles.

The correct strategy would be to apply the Drift operation to all particles, then the Kick operation to all the particles, then the Drift operation to all the particles again.

The difference is the same as the difference between the following two statements:

  • For each particle i, for each integration step k apply step k to particle i.
  • For each integration step k, for each particle i apply step k to particle i.

Now it looks like someone interchanged the order of two nested for loops. Unlike matrix multiplication where reordering of the loops doesn’t change the end result, for an n-body integrator (like our particle system) the reordering of the loops has a significant impact on the correctness (and stability) of the end result.

This arises due to the fact that matrix multiplication doesn’t have dependencies on the other entries in the matrices whilst in a particle system there are interdependencies between the particles (especially in the case of SPH/Fluid Simulations).

I may be wrong about this, and so I’m posting it here first before prodding a developer and seeking answers.

(NOTE: RK4 should also see the same issue as it has 4 substeps. Since it isn’t energy conservative you shouldn’t be using it for fluid or newtonian particle systems, it can lead to some very odd behaviour like particles randomly exploding)