Two Falling Dots

The galaxy at the top of this page is a live simulation: 40,000 bodies, every one pulling on every other, stepped forward every frame on your graphics card. It runs in slow motion so you can watch it turn. The last demo on this page runs the same code with up to 300,000.

Done the obvious way, 300,000 bodies means 300,000² pair forces per step. That's 90 billion, and in a browser tab it would take more than a minute to draw a single frame. This article is how I got it running in real time. Each step fixes whatever was limiting the step before it, and each one has a demo you can play with.

Step 1 · One forceThe formula

All of the physics in this article is one line. Everything after this section is about evaluating it faster.

// acceleration of body i, caused by body j
d   = position[j] - position[i]
a_i += G * mass[j] * d / (|d|² + ε²)^(3/2)

That's Newton's law of gravity with one fudge added, the ε² term, called softening. Real gravity goes to infinity as two points touch. In a simulation that moves in fixed time steps it's worse than that: a body passes close to another, the force spikes for a single step, and the body gets thrown out of the system faster than anything real could move. Softening caps the force. Inside a distance of about ε it stops growing, as if each point mass had been smeared into a small ball of radius ε.

Try it below. Press on the black, drag in the direction you want the body to go, and let go. Then pull the softening slider all the way down and launch a body straight at the star. It comes back out far faster than it went in. That spike is what ε is there to prevent.

One more detail matters here, and it costs nothing. Update the velocity from the acceleration first, then move the body using the new velocity. This is called semi-implicit Euler. Do it the other way around, moving the body with its old velocity first, and orbits slowly gain energy and spiral outward. The demo has a button for this. Watch the inner planet: switch to "position first" and it starts to unwind within a few laps.

One star and a few planets is cheap, because each planet only feels the star. It gets expensive when every body feels every other body.

Step 2 · All pairsEveryone pulls on everyone

A galaxy has no central star holding it together. The disk holds itself together with its own gravity, so every body pulls on every other one. With three or more bodies there's no formula for what happens next. The only way to find out what a thousand gravitating bodies will do is to step them forward and watch.

With a computer it's two nested loops:

for (let i = 0; i < count; i++) {
  let ax = 0, ay = 0;
  for (let j = 0; j < count; j++) {
    const dx = x[j] - x[i], dy = y[j] - y[i];
    const r2 = dx*dx + dy*dy + eps2;
    const inv = G * mass[j] / (r2 * Math.sqrt(r2));
    ax += dx * inv;  ay += dy * inv;
  }
  vx[i] += ax * dt;  vy[i] += ay * dt;
}

This is correct, and it's the reference everything later gets checked against. It's also slow. Drag the slider right and watch the milliseconds per step: the work grows with the square of the body count. Twice the bodies is four times the work. 4,000 bodies is 16 million force calculations per step, at sixty steps a second.

Somewhere along that slider your machine ran out of room. The results are still right, and the inner loop is about as tight as JavaScript gets. An O(n²) loop outgrows any constant speedup eventually. To go further you need a lot more arithmetic or a better algorithm, and I ended up needing both. The extra arithmetic is already sitting in your machine.

Step 3 · ParallelMove it to the GPU

Look at the nested loop again. Body 0's force doesn't depend on body 1's, so every pass of the outer loop is independent. That's the kind of problem a GPU is built for: thousands of identical calculations that never need to talk to each other.

So the GPU version runs one thread per body, and each thread does the inner loop. The outer loop goes away. The hardware runs all of its iterations at once.

The most important thing about GPUs shows up right away: they usually run out of memory bandwidth before they run out of math. If 100,000 threads each read all 100,000 positions from GPU memory, the memory system can't keep up. So the threads work in groups of 256. A group loads a tile of 256 positions into fast on-chip shared memory, every thread sums forces against that tile, and then the group loads the next one.

var<workgroup> tile: array<vec3f, 256>;  // fast on-chip memory

for (var t = 0u; t < numTiles; t++) {
  tile[local_id] = load(bodies[t * 256u + local_id]); // each thread fetches one
  workgroupBarrier();                                  // wait for the whole tile
  for (var k = 0u; k < 256u; k++) {
    acc += force(tile[k], myPosition);                 // everyone uses all 256
  }
  workgroupBarrier();
}

Each position now comes out of main memory once per group instead of once per thread, which cuts memory traffic by 256×. This is the tiled n-body kernel from GPU Gems 3, and it's more or less the hello world of GPU computing.

Same slider as before, and the wall that was around 4,000 bodies is now dozens of times further out. The problem underneath hasn't changed, though. Watch the readout: at 120,000 bodies you're computing 14.4 billion pair forces per step, and that still grows with the square of the count. The GPU moved the wall. Getting past it means computing fewer forces, not the same forces faster.

Step 4 · ApproximationGroup the far-away bodies

Every star in the Andromeda galaxy is pulling on you right now, and you feel all of them as one pull toward one smudge in the sky. From far enough away, a cluster of bodies acts like a single point with the cluster's total mass, sitting at its centre of mass.

Josh Barnes and Piet Hut turned that into an algorithm in 1986. Build a quadtree: split space into four squares, split any square holding more than one body into four again, and keep going. Each cell stores its total mass and centre of mass. To get the force on a body, start at the root and ask one question at every cell:

// is this cell small enough, seen from here, to count as one point?
if (cellWidth² < θ² × distance²) {
  acc += forceFrom(cell.totalMass, cell.centerOfMass);  // one term for the whole cell
} else {
  descend into the four children;                       // too close, look closer
}

θ (theta) trades accuracy for speed. It sets how small a cell has to look, its width divided by its distance, before it counts as a single point. Cells near the body get opened all the way down to individual bodies, and the far side of the galaxy collapses into a few lumps. A sum that took n terms now takes about log n.

Here's one walk, a piece at a time. Scroll slowly, and scroll back up whenever you want to see a stage again.

And here it is live. Your cursor is the test body, and each yellow box is a cluster the walk treated as one point mass for that spot. Move toward the middle of the disk and the boxes near you shrink while the far ones stay big. Then drag θ down to 0.1 (lots of small boxes) and up to 2.0 (a handful of boxes for the whole galaxy).

The readout checks the walk against the exact sum. At the default θ = 0.8 it uses about 60 terms instead of 3,000 and usually lands within a few percent of the exact force. At θ = 0.3 the error drops below 1%, for a few hundred terms. That's still a tenth of the work.

The obvious next move is to port the tree walk to the GPU and add a zero to the body count. One thing gets in the way, and it's the tree.

Step 5 · The bottleneckA tree with nothing to build

The textbook Barnes–Hut tree is built one body at a time. Find the leaf the body belongs in, and if something's already there, split the leaf and push the old occupant down a level. Each insert depends on the tree the last one left behind, so the build is sequential. And because every body moves every step, the tree has to be rebuilt from scratch every step.

On a CPU that's fine at small counts. But the simulation lives on the GPU now, so every frame would mean copying all the positions back to the CPU, building the tree in single-threaded JavaScript, and uploading it again. Past about 100,000 bodies, that round trip costs more than the physics. There's nothing wrong with the algorithm. The tree is being built on the wrong side of the bus.

You can port the sequential build to the GPU with atomics and locks. People have, and it gets complicated. I went the other way and changed the tree so there's nothing left to build.

Make the quadtree complete: split every cell to the same depth, whether anything is in it or not. Then the whole structure is arithmetic. The children of cell (level, ix, iy) are always (level+1, 2ix+dx, 2iy+dy), so there are no pointers to allocate or follow, and each level is a flat block in one buffer. If you've worked with texture mipmaps, you already know this structure. It's a mipmap of mass.

Slide through the levels below. Level 7 is a 128×128 grid over the live disk, and each level above it merges 2×2 cells into one, up to a single cell holding the whole galaxy. Each cell stores what a Barnes–Hut node stores, total mass and centre of mass (the pink dots), so the walk from the last demo runs on it unchanged.

Unlike the pointer tree, this one builds in parallel. It takes two passes.

Scatter. Every body adds its mass into its cell at the finest level, all at the same time. Thousands of threads hit the same cells, so the adds have to be atomic:

// one thread per body, in any order, nobody waits
atomicAdd(&grid[cell],     u32(round(m * FP_SCALE)));                     // Σ mass
atomicAdd(&grid[cell + 1], bitcast<u32>(i32(round(m * dx * FP_SCALE)))); // Σ m·dx
atomicAdd(&grid[cell + 2], bitcast<u32>(i32(round(m * dy * FP_SCALE)))); // Σ m·dy

Reduce. Build each coarser level from the one below it: every parent adds up its four children. One thread per parent, no conflicts, one pass per level. A 1024×1024 pyramid takes ten of these passes and about a millisecond in total. The sequential tree build has turned into one scatter and a mipmap chain.

Step 6 · GPU-residentThe whole pipeline on the GPU

Put it together and each simulation step is a fixed list of GPU dispatches, fifteen at full size, and none of them touch the CPU:

clear    → zero the grid's atomic accumulators
bounds   → where is everything? (parallel min/max)
scatter  → every body adds its mass to its finest cell   (the atomics above)
resolve  → fixed point back to mass + centre of mass
reduce   → build the pyramid, one level per pass         (10 passes at 1024²)
force    → every body walks the implicit tree            (Barnes–Hut, no pointers)
           …and moves itself, in place

Two details are worth pointing out. The bounds pass refits the grid to the bodies every step, so the pyramid always covers everything, even when the galaxy drifts or flings bodies out. And the force pass writes each body back in place. That's safe because no kernel ever reads another body's position directly; everything a body knows about the others comes through the pyramid. So there's one body buffer and no second copy to swap with.

Here's the result. It starts at 100,000 bodies, and the θ slider is the same one from the 3,000-body demo. Steps per frame runs more simulation between frames, which speeds up time without asking more of your display.

On my machine this runs 300,000 bodies at 60 frames a second with two steps per frame, about 36 million body-steps a second. The nested loop from Step 2 would need over a minute for one step at that size. The physics hasn't changed: same formula, same softening, and a θ you can check against the smaller demos. Here's where the speedup came from:

two loops in JS             ~4,000 bodies          (Step 2)
+ GPU, one thread per body   dozens of times more   (Step 3)
+ tiles in shared memory     (part of that)         (Step 3)
+ Barnes–Hut θ test          O(n²) → O(n log n)     (Step 4)
+ pointer-free pyramid       tree build ~1 ms       (Step 5)
+ nothing leaves the GPU     no round trips         (Step 6)

This demo and the galaxy at the top have one thing the earlier ones don't: a dark matter halo. It's a single fixed pull toward the centre, added to each body's acceleration after the tree walk. Turn it off with the dark halo button and the disk can't hold itself together. It breaks into clumps, the clumps settle into a ring, and bodies get thrown clear. Turn it back on and the disk stays put and grows spiral arms instead. Colour is still speed: the slow-turning core shows up gold and the faster disk blue-white.

EpilogueWhat carries over

The gravity is specific to this page, but a few lessons run through the whole series:

• Figure out which wall you're hitting. The GPU moved the O(n²) wall. Only Barnes–Hut got rid of it.
• Accuracy is a dial. θ gives up a few percent of accuracy for a huge cut in work, and you get to choose where it sits.
• Data layout decides what can run in parallel. The pointer quadtree and the pyramid hold the same information. One can't be built in parallel and the other is a scatter and a mipmap.
• Watch the round trips. The CPU tree build didn't lose because CPUs are slow. It lost because of the trip across the bus.

Things I'd try next: sort the bodies along a space-filling curve (Morton order) so neighbouring GPU threads walk the same parts of the tree; replace the uniform pyramid with an adaptive tree built in parallel (an LBVH) for very clumpy scenes; and add a third dimension. Part two does that last one.

The galaxy at the top of the page is still running. Scroll back up and you can name everything it's doing.

Next: Part two, Gravity on a Dome. The simulation goes 3D. Quadtrees become octrees, and the galaxy gets bent over a hemispheric dome.