Grouping the particles

Up a level : Themodynamics
Previous page : Colliding non-rotating balls in 2D
Next page : Motion in the ocean… I mean in the atmosphere

If we have many particles and if use the algorithm from the previous page we will sooner or later notice that the calculations will take too long to make a smooth simulation. The reason for this is that we need to check every pair of particles. If we have for example 3000 particles, we will need to do 2999+2998+2997+…+3+2+1 tests for collisions between particles.

We can find the sum by adding the terms of this series with the same series in revered order. This will give us 2999 terms with the value 3000 since 2999+1=3000, 2998+2=3000 and so on.

This will thus give us 2999·3000/2=4498500 comparisons.

Now let us instead split the area of the simulation in say vertical strips, and that each strip is as wide as the diameter of the particle.

If two particles can collide they must be in the same strip or the adjacent strips.

We can then loop through all particles, and place the identity of each particle that has its centre in a slot in a bin (in this case an array) for each strip. That would be 3000 comparisons. The good thing now is that we need only to check each particle with the particles in its own bin and the two neighbouring bins. In practice we test all particles in one bin with all others in the same bin, then each particle with all bins to the right. Then we do the same things with the next bin to the right, and so on until the last bin, where we of course do not need to check it against the next to the right, since we have no more bins.

Say we have balls with a diameter of 4 pixels, and the total width of the area is 600 pixels, then we need 600/4=150 bins. On average we will get 3000/150=20 particles per bin. That would give us, on average, 20·19/2 compassions per bin, or 190·150=28500 compassions. Next, we need to check against all particles in the bin to the right. That would give us 20·20·149= 59600 compassions. This adds up to 91100 compassions. That would be a factor of over 49 times fewer comparisons.

The time required for the actual collisions will remain the same since we will have the same number of collisions. The difference is in the number of comparisons to find out if it even is a collision.

So how much will this change the times in practice? On the laptop I used, the time for the naïve approach was 36 ms for 3000 particles – and that includes collisions with the walls, with each other, and drawing the particles. Since the computer I used has a framerate of 60 frames per second, i.e. about 17 ms per frame, we will use two frames per calculation, and our frame rate will drop to 30 frames a second, which is noticeable indeed. With the above algorithm, the time falls to 11.5 ms, and we will be able to do one full round of calculations within one frame. The time for the collision checks and calculations falls from 28 ms to 3.4 ms.

An alternative arrangement

How about arranging the slots in the other direction? On one hand some lower slots contain way more particles than before (one might guess), but on the other hand, the upper ones contain basically at most one each – and also, the number of slots will, in this case, be way more since the simulation area is 6000 pixels high. OK, so let’s test this.

We get a total time of about 10.4 ms, and the time for the collision calculations goes down to about 2.0 ms, so this may seem to be the way to go – even though the difference was not that big.

It is harder to calculate the expected number of comparisons now since the number of particles per slot will vary widely from the average, and since we cannot use the average, since the number of comparisons basically depends on the number of particles squared. What I did though was to add a bit of code that counted the number of comparisons, and it seems to end up at about 70000 to 80000.

The total running time for one turn of the average varies more with the horizontal slots, so I will use the vertical slots though.

Here is the new collision function

// This is to check if particle i and j
// has collided. The horizontal distance between
// the particles is checked to see if a collision
// is possible, if so, then the actual distance is
// checked. If a collision is still
// possible then it is checked if the particles are moving
// toward each other => collision, or not.
// Then the calculations to simulate
// a totally elastic collision is done.
// This is free to use - but I'm happy if you add a link to 
// this page in your code. 
function collision(i, j) {
  let dy = particles[i].y - particles[j].y;
  if (Math.abs(dy) < diameter) {
    let dx = particles[i].x - particles[j].x;
    let dy = particles[i].y - particles[j].y;
    let d2 = dx ** 2 + dy ** 2;
    if (d2 < diameter2) {
      let dvx = particles[j].vx - particles[i].vx;
      let dvy = particles[j].vy - particles[i].vy;
      let dvs = dx * dvx + dy * dvy;
      if (dvs > 0) {
        dvs = dvs / d2;
        dvx = dvs * dx;
        dvy = dvs * dy;
        particles[i].vx += dvx;
        particles[i].vy += dvy;
        particles[j].vx -= dvx;
        particles[j].vy -= dvy;
      }
    }
  }
}

// This is to check if particle i and j

// has collided. The horizontal distance between

// the particles is checked to see if a collision

// is possible, if so, then the actual distance is

// checked. If a collision is still

// possible then it is checked if the particles are moving

// toward each other => collision, or not.

// Then the calculations to simulate

// a totally elastic collision is done.

// This is free to use - but I'm happy if you add a link to

// this page in your code.

function collision(i, j) {

let dy = particles[i].y - particles[j].y;

if (Math.abs(dy) < diameter) {

let dx = particles[i].x - particles[j].x;

let dy = particles[i].y - particles[j].y;

let d2 = dx ** 2 + dy ** 2;

if (d2 < diameter2) {

let dvx = particles[j].vx - particles[i].vx;

let dvy = particles[j].vy - particles[i].vy;

let dvs = dx * dvx + dy * dvy;

if (dvs > 0) {

dvs = dvs / d2;

dvx = dvs * dx;

dvy = dvs * dy;

particles[i].vx += dvx;

particles[i].vy += dvy;

particles[j].vx -= dvx;

particles[j].vy -= dvy;

}

This is the particle-particle interaction checker:

// In this function particle-particle interactions
// is checked. First the particles will
// be placed in vertical slots with a width
// equal to (or slightly larger than) the
// radius of the particles.
// Then we loop over the slots, checking for
// collisions among particles in the same slot,
// and then with particles in the slot to the left.
// I.e. the slots are checked all in 0 against all in 0,
// all in 0 against all in 1, all in 1 against all in 1,
// i.e. 0-0, 0-1, 1-1, 1-2, 2-2 and so on until the last slot
// where we of course do not check against the non existent
// next slot.
// In this way we ensure that all possible pairs of particles
// that could possibly collide are checked.
// This is free to use - but I'm happy if you add a link to 
// this page in your code. 
function particleParticle() {
  slots = Array(nSlots)
    .fill(null)
    .map(() => []);
  // Place the particles in the correct slot.
  for (let i = 0; i < n; i++) {
    let s = Math.floor(particles[i].x / diameter);
    s = Math.max(s, 0);
    s = Math.min(s, nSlots - 1);
    slots[s].push(i);
  }
  for (let s = 0; s < nSlots; s++) {
    // Check all pairs of particles in the current slot.
    for (let i = 0; i < slots[s].length - 1; i++) {
      for (let j = i + 1; j < slots[s].length; j++) {
        collision(slots[s][i], slots[s][j]);
      }
    }
    // Check all particles in the current slot
    // against all particles in the slot to the right
    // unless we are at the last slot.
    if (s !== nSlots - 1) {
      for (let i = 0; i < slots[s].length; i++) {
        for (let j = 0; j < slots[s + 1].length; j++) {
          collision(slots[s][i], slots[s + 1][j]);
        }
      }
    }
  }
}

// In this function particle-particle interactions

// is checked. First the particles will

// be placed in vertical slots with a width

// equal to (or slightly larger than) the

// radius of the particles.

// Then we loop over the slots, checking for

// collisions among particles in the same slot,

// and then with particles in the slot to the left.

// I.e. the slots are checked all in 0 against all in 0,

// all in 0 against all in 1, all in 1 against all in 1,

// i.e. 0-0, 0-1, 1-1, 1-2, 2-2 and so on until the last slot

// where we of course do not check against the non existent

// next slot.

// In this way we ensure that all possible pairs of particles

// that could possibly collide are checked.

// This is free to use - but I'm happy if you add a link to

// this page in your code.

function particleParticle() {

slots = Array(nSlots)

.fill(null)

.map(() => []);

// Place the particles in the correct slot.

for (let i = 0; i < n; i++) {

let s = Math.floor(particles[i].x / diameter);

s = Math.max(s, 0);

s = Math.min(s, nSlots - 1);

slots[s].push(i);

}

for (let s = 0; s < nSlots; s++) {

// Check all pairs of particles in the current slot.

for (let i = 0; i < slots[s].length - 1; i++) {

for (let j = i + 1; j < slots[s].length; j++) {

collision(slots[s][i], slots[s][j]);

}

// Check all particles in the current slot

// against all particles in the slot to the right

// unless we are at the last slot.

if (s !== nSlots - 1) {

for (let i = 0; i < slots[s].length; i++) {

for (let j = 0; j < slots[s + 1].length; j++) {

collision(slots[s][i], slots[s + 1][j]);

}

Up a level : Themodynamics
Previous page : Colliding non-rotating balls in 2D
Next page : Motion in the ocean… I mean in the atmosphere

Last modified: Dec 28, 2023 @ 18:03