Gravity on a Dome — Algorythmic Explorations

Part one ended with 300,000 bodies on a flat plane, using a Barnes–Hut tree with no pointers that the GPU rebuilds from scratch every step. This part takes the same solver into 3D, and then restricts the 3D world to the surface of a dome.

The dome at the top of the page is 20,000 bodies on a thin hemispheric shell, stepped twice per frame. They pull on each other through 3D space while a spring keeps them on the surface. I halved gravity for that one so it keeps flowing instead of settling into clumps. Drag to orbit it, and pinch or hold ⌘/Ctrl while scrolling to zoom (plain scrolling scrolls the page).

The 3D port turned out to be small. Memory is the part that changes most, and the dome is one extra force term.

Step 1 · Into 3DQuadtree to octree

Going from a quadtree to an octree changes very little. The only structural difference is how a cell finds its children:

// 2D: children of (level, ix, iy), four of them
(level+1, 2ix+dx, 2iy+dy)                  dx, dy ∈ {0, 1}

// 3D: children of (level, ix, iy, iz), eight of them
(level+1, 2ix+dx, 2iy+dy, 2iz+dz)          dx, dy, dz ∈ {0, 1}

Children are still computed rather than stored, and each level is still a flat block in one buffer. The Barnes–Hut test, cellWidth² < θ² × distance², doesn't change at all, since widths and distances are lengths in any number of dimensions.

What does change is memory. Each level now grows as a cube instead of a square. Part one's finest grid went up to 1024², about a million cells and 16 MB of accumulators. The same depth in 3D would be 1024³, a billion cells and around 16 GB, which isn't going to happen in a browser tab. So I capped the octree at seven levels: 128³, about two million cells and 34 MB. Below that cap the solver picks its depth from the body count, ceil(log₂(n) / 3) clamped to between 4 and 7, which gives roughly as many finest cells as bodies.

The demos on this page never reach the cap. At the default 20,000 bodies the finest grid is 32³, and at the slider's maximum of 200,000 it's 64³. It would take more than 2¹⁸ = 262,144 bodies to get a 128³ grid.

A coarser floor costs resolution up close. The finest level is where the walk stops, so the other bodies in each nearby finest cell get lumped into one point at that cell's centre of mass. At 20,000 bodies the finest cells start out about 0.06 wide, close to the softening length of 0.05, so the grid blurs about as much as softening already does. For a galaxy, where no single close encounter matters, you can't see the difference. For planets making close passes it would be the first thing I'd fix.

The bookkeeping that changes from 2D to 3D

Level offsets. The levels above level l hold 1 + 8 + 64 + … + 8^l−1 nodes, which sums to (8^l−1)/7. In 2D it was (4^l−1)/3. Both denominators are 2^d−1. It also means all the levels above the finest one together add only about a seventh more memory in 3D (a third in 2D), so once the bottom level is built the rest of the tree is cheap.

The spare word gets used. Part one's scatter used a stride of four atomic words per cell, because power-of-two strides keep the indexing cheap, but only filled three: mass, Σm·dx and Σm·dy. In 3D the fourth word holds Σm·dz, and the layout is full:

atomicAdd(&grid[c],     u32(round(m * FP_SCALE)));  // Σ mass
atomicAdd(&grid[c + 1], …(m * frac.x * FP_SCALE));  // Σ m·dx
atomicAdd(&grid[c + 2], …(m * frac.y * FP_SCALE));  // Σ m·dy
atomicAdd(&grid[c + 3], …(m * frac.z * FP_SCALE));  // Σ m·dz, the word 2D left empty

Bounds and stack. The bounds pass fills six slots (min and max of x, y and z) instead of four. Stack entries pack into 32 bits as level(4) | iz(9) | iy(9) | ix(9). Opening a cell now pushes eight children instead of four, so the stack has to be deeper. The shader uses a fixed 64-entry array, enough for the seven-level cap. At the default 32³ grid the walk never needs more than about 36 entries.

Mass goes in position.w. A 3D position needs three floats and GPU buffers want 16-byte alignment, so there was going to be a fourth float anyway. Putting the mass there removes the separate mass buffer. A body is one vec4 of position and one of velocity, 32 bytes and one fetch each, so 200,000 bodies fit in about 6 MB.

Step 2 · RenderingDrawing without sorting

Before the octree is any use I need to be able to see in 3D, and drawing transparent things in 3D usually means sorting them. Normal blending depends on order: draw a near glow before a far one and the picture comes out wrong. Sorting a few hundred thousand particles by depth every frame would be a whole extra pipeline.

I skip the sort by having the particles add light instead of blocking it. Addition doesn't care about order (a + b = b + a), so any number of overlapping glows add up to the same picture whichever order they're drawn in. There's no sort and no depth buffer. Stars mostly work this way anyway. They give off light, and at galactic distances they almost never cover each other.

Aside · Fifty years of draw order

Draw order has been a problem for about as long as there's been computer graphics. The painter's algorithm from the early 1970s does what the name says: draw the far things first and paint the near ones over them. The z-buffer, described by Ed Catmull and independently by Wolfgang Straßer in 1974, removed the sort for opaque scenes by storing a depth for every pixel, and every GPU now has one in hardware. A depth test can only keep or discard a pixel, though, so it doesn't help with blending, and transparency brought the sort back. Order-independent transparency is still an active research topic. A field of stars is one of the few scenes where plain additive blending is a fair model.

The geometry is cheap too. Each body is drawn as a small square facing the screen, and there's no vertex buffer. The six corners are hard-coded in the vertex shader, which looks up each body's position and velocity by instance index, straight from the buffers the simulation writes. Nothing is copied between the physics and the pixels. The fragment shader fades each square out from the centre, and speed picks the colour: blue for slow bodies, orange for fast ones. Bodies in a collapsing clump move fast, so clumps glow orange. The light adds up in a 16-bit float buffer, so dense regions keep getting brighter instead of clipping, and a final pass tone-maps that onto the screen.

Below is part one's tiled O(n²) kernel with a z component added. Drag to orbit it. The disk has some thickness now, so you get parallax.

It starts at 25,000 bodies and the slider tops out at 60,000. Brute force costs n² in any number of dimensions, which is what the octree is for.

Step 3 · The octreeRebuilding the octree every step

The 3D pipeline is part one's pipeline with the constants changed. Every body moves every step, so every step rebuilds the octree from nothing. At the default size that's ten dispatches:

clear     → zero the grid's atomic words, reset the bounds
bounds    → min/max over every body                 (the root cell)
scatter   → each body adds into its cell's 4 words  (fixed-point atomics)
resolve   → fixed point → centre of mass and mass   (finest level)
reduce ×5 → each level from its 8 children          (×6 at 64³)
force     → stack walk, θ test, integrate           (in place)

At 20,000 bodies the grid is 32³, so the clear pass zeroes about 131,000 words and there are five reduce passes. At 200,000 bodies it's 64³, about a million words and six passes. As in 2D, the force kernel reads only its own body from the body buffers (everything it knows about the others comes from the octree), so it integrates in place, and the renderer draws from that same buffer.

How a body avoids pulling on itself

Every body scatters its own mass into the grid, so every cell on the path from the root down to the body's finest cell includes the body itself. If the walk accepted one of those cells as a single point, the body would be pulled toward a centre of mass that is partly itself. That can happen: a big cell can pass the θ test when most of its mass sits at the far end of the cell from the body.

So the walk has two rules. At the coarse levels it never accepts the cell the body sits in. It always opens that one, all the way down. At the finest level, when it reaches the body's own cell, it subtracts the body's mass from the cell's total and takes the body's position back out of the centre of mass, then attracts toward whatever is left. If the body was alone in the cell, nothing is left and the cell is skipped.

// coarse levels: never summarize the cell this body sits in
let shift = FINEST - lvl;
let mine = ix == (myIx >> shift) && iy == (myIy >> shift) && iz == (myIz >> shift);
if (!mine && w * w < theta2 * r2) { /* accept as one point */ } else { /* push 8 children */ }

// finest level, own cell: take this body back out
m = n.w - myMass;
if (m <= 1e-9) { continue; }             // alone in the cell
com = (n.xyz * n.w - p * myMass) / m;

It starts at 20,000 bodies, and the slider goes from 10,000 to 200,000. θ runs from 0.4 to 1.4 and starts at 0.8, the same default as part one. The starting galaxy is part one's recipe with some thickness: a spinning disk where each body gets the circular velocity for the mass inside its radius, plus a little vertical scatter. Leave it running and it does what the flat disk did in part one. It breaks into clumps, and the clumps settle into a ring. That's the same cold-disk instability from part one's aside, and none of it is scripted.

Step 4 · ConstraintsThe dome

Everything so far is gravity in open space, where bodies go wherever the forces send them. For the dome I wanted the bodies to live on a surface instead, a hemisphere they flow across like weather over the top of a planet. You could do that with spherical coordinates, or with an integrator written for curved surfaces. I went with a spring, which is a few lines in the force kernel:

if (P.shellK > 0.0) {
  let r = length(p);
  let rhat = p / max(r, 1e-6);
  acc = acc - rhat * (r - P.shellR) * P.shellK;  // radial spring: stay on the shell
  if (p.z < 0.0) {
    acc.z = acc.z - p.z * P.shellK * 4.0;        // floor spring: stay above the equator
  }
}

A body inside the shell radius gets pushed out, one outside gets pulled in, and one that drops below the equator gets pushed back up. I made the floor spring four times stiffer than the shell spring, because I wanted the equator to act like a hard edge and the shell like a soft preference. Each push is proportional to how far the body is from where it should be. This is called a penalty force. The constraint never edits positions and never touches the integrator or the octree. It's one more term in the acceleration, next to gravity.

Below is one body on the spring, in a 2D slice through the dome. It's a real simulated trajectory, computed once when the page loads, and your scroll position scrubs back and forth through it.

The demo below starts with the dome on: 20,000 bodies spread over the hemisphere, dome strength 6, and a little velocity damping (0.999 per step) so the shell settles instead of ringing. Drag the strength down to zero and the springs let go. Within a few seconds free 3D gravity pulls the whole shell together into one tight clump. Then raise the strength again. The springs drag the clump back out to the shell, but it stays one lump and slides around on the surface; it doesn't spread back out into a dome. Press re-seed to get the shell back.

While the bodies are held on the surface, gravity still acts through 3D space. Two clumps on opposite sides of the dome pull on each other along the straight line through the inside of the dome, not along the surface. The spring cancels the part of that pull that points into the dome, and what's left slides each clump along the shell toward the other.

Why a spring and not a hard projection

The obvious alternative is projection: after each step, set p = normalize(p) * shellR. Projection moves the body without a force, so the integrator never sees it. You either throw away the radial velocity by hand as well, or leave it in and have the body fight the projection every step. Either way energy changes without anything in the simulation accounting for it. A penalty force stays inside the integrator, where damping and gravity act on it like any other force and I can tune it. The small velocity damping that switches on with the dome is there so the shell can settle instead of oscillating forever.

The downside of a spring is stiffness. An explicit integrator can only handle a spring so stiff before the spring's own oscillation is faster than the time step can follow, and the simulation blows up. That's why the dome uses a modest strength plus damping rather than a rigid wall. If penalty stiffness ever became the limit, the next things to try would be Lagrange multipliers or the position-based dynamics that a lot of game physics uses. For this, a spring is enough.

Nothing in the spring cares that the shape is a sphere, either. length(p) − shellR is the signed distance to the shell, and the spring pushes along its gradient. Swap in a different signed distance function, like a torus or a terrain heightfield, and gravity would flow over that shape instead, held there by the same few lines.

EpilogueWhat changes in 3D

A few things from this part that I'll keep using:

An extra dimension mostly costs memory. Going from 2D to 3D changed constants, 4 to 8 and squares to cubes, but none of the algorithms. The one real decision was how deep a grid I could afford, and the answer was a coarser floor that softening mostly covers.
Operations that ignore order save work. Additive blending gives the same result in any draw order, so there's no depth sort to write.
Constraints can be forces. The dome is a few lines in the force kernel, and nothing else in the solver knows it's there. Part three uses the same idea for the walls of a water tank.
Physics and rendering can share a buffer. The renderer reads the simulation's buffers directly, by instance index, in the same frame they were written.

Things I'd try next: sort the bodies in Morton order so neighbouring GPU threads walk the same branches of the octree; replace the uniform pyramid with an adaptive tree (an LBVH) so very clumpy scenes don't run out of resolution at the finest level; and add collisions, which would turn the dome into a container for particles that bounce off each other instead of passing through.

Next: Part three, Sorting Water. The particles stop attracting and start colliding, and finding each one's neighbours fast turns into a sorting problem.