01 — motivation
What "GPU compute" even means, and why a fluid needs it
A CPU is a small team of very fast, very clever workers. It has a handful of cores, and each one races through a list of instructions one after another. If you ask a CPU to update 20,000 water particles, it walks the list: particle 0, then 1, then 2, all the way to 19,999. Fast, but fundamentally one thing at a time per core.
A GPU is the opposite shape: thousands of much simpler workers that all do the same instructions at the same time, each on a different piece of data. That is exactly the shape of a particle simulation — every particle runs the identical update rule, just with its own position and velocity. So instead of a loop that visits 20,000 particles in sequence, you launch 20,000 tiny threads and they all move in parallel. This is called a compute shader: a program that doesn't draw anything directly, it just crunches numbers in parallel and leaves the results in memory.
The fluid in this project is an SPH simulation — Smoothed Particle Hydrodynamics. The water isn't a surface; it's 20,000 little balls that push on their neighbours. Every frame, each particle has to look at the particles near it, feel how crowded it is, and nudge itself to relieve the pressure. Done on a CPU, that's 20,000 particles × dozens of neighbours × several times per frame — far too slow for a phone. Done on a GPU, all 20,000 are handled together. That is why this lives on the GPU.
02 — the core idea
Write once, run on both
The heart of zimr's compute system is a module called kompute. You write your simulation logic one time, in plain Zig. From that single source, the build produces two things:
- A native CPU version — an ordinary loop you can run, step through in a debugger, and trust.
- A GPU version — the same logic compiled all the way down to a WebGPU compute shader.
How can one function become both? Zig has a value, available at compile time, that says which machine it's being compiled for. kompute exposes it as is_gpu:
/// True when compiling for SPIR-V (the GPU dispatch). False for the
/// native CPU loop. Kernel files comptime-branch on it.
pub const is_gpu: bool = builtin.target.cpu.arch.isSpirV();Because is_gpu is known at compile time, an if (is_gpu) isn't a runtime cost — the compiler keeps one branch and deletes the other. The CPU build and the GPU build literally contain different code, generated from the same text. The CPU version becomes the oracle: it is the definition of "correct." If the GPU version ever disagrees with it, the GPU version is wrong. That single idea — a trustworthy CPU twin you can compare against — is what makes it possible to develop a parallel program without going mad.
The same compile-time switch is how kompute handles operations the GPU compiler can't express directly. Take atomics — the special "no two threads collide" operations the particle grid needs (more on those soon). The CPU has real atomic instructions; the GPU shader language has its own. So the helper picks the right one at compile time:
/// Atomic read-modify-write add. Returns the value stored BEFORE the
/// add (WGSL atomicAdd semantics — the grid build uses the returned
/// value as the claimed slot index).
pub inline fn atomicAdd(arr: anytype, idx: u32, val: u32) u32 {
if (is_gpu) {
return zatomicAdd(arr, idx, val); // GPU: rewritten into WGSL later
}
return @atomicRmw(u32, &arr[idx], .Add, val, .monotonic); // CPU: the real thing
}One call, k.atomicAdd(...), written once in the simulation. On the CPU it becomes a genuine atomic instruction; on the GPU it becomes a placeholder that a later stage turns into the WGSL atomic builtin. The kernel author never sees the difference.
03 — anatomy
What a kernel file actually contains
A "kernel" is one parallel function. A kernel file declares three things, and then one or more kernel functions. Here's the shape, straight from kompute's own documentation:
const k = @import("kompute");
// 1. config — how big, and how many threads per workgroup
pub const config = k.Config{ .max = 1024, .workgroup = 64 };
// 2. Buffers — the data living on the GPU (one big struct)
pub const Buffers = extern struct { data: [config.max]f32 };
// 3. Params — small read-only values passed in each dispatch
pub const Params = extern struct { count: u32, _pad: [3]u32 = .{0,0,0} };
pub const g = k.Globals(@This()); // wires up storage on both targets
const b_data = g.bind(.data); // one alias per Buffers field
// the kernel: runs once per thread, c.id is "which thread am I"
pub fn double(c: k.Ctx(@This())) void {
if (c.id >= c.params.count) return; // ignore threads past the data
b_data[c.id] = b_data[c.id] * 2.0;
}
comptime { k.installKernel(@This(), "double"); } // register itThree pieces deserve a closer look.
Buffers — the shared data
Buffers is one big struct holding everything the simulation keeps on the GPU: arrays of positions, velocities, and so on. It uses extern struct so its memory layout is fixed and predictable — the CPU and the GPU have to agree, byte for byte, on where each field sits.
Each field is reached through an alias made by g.bind(.field). On the GPU these become separate storage bindings; on the CPU they're pointers into plain memory. Either way, b_data[i] means the same thing in both worlds.
Ctx — what each thread knows
When a kernel runs, it receives a tiny context: c.id (this thread's global index — "which particle am I?") and c.params (the read-only settings for this dispatch). That's deliberately almost everything a thread needs and nothing it doesn't:
pub fn Ctx(comptime Module: type) type {
return struct {
id: u32, // this invocation's global index
params: Module.Params, // the kernel's Params, by value (read-only)
};
}The guard if (c.id >= c.params.count) return; appears at the top of every kernel. The GPU launches threads in fixed-size groups, so you almost always get a few more threads than you have data. Those extra threads must do nothing. Forgetting this guard is how you read or write past the end of an array on the GPU.
installKernel — making it launchable
The comptime { k.installKernel(...) } line is what turns your function into something the GPU can actually start. On the GPU build it generates the official entry point: read this thread's index, build the Ctx, call your function. On the CPU build it does nothing (the host just calls your function in a loop). Here it is:
pub fn installKernel(comptime Module: type, comptime name: []const u8) void {
if (!is_gpu) return; // CPU: nothing to generate
const Entry = struct {
fn run() callconv(.spirv_kernel) void {
@setRuntimeSafety(false);
const c: Ctx(Module) = .{
.id = gpu.global_invocation_id[0], // ask the GPU "who am I"
.params = Module.g.P,
};
@field(Module, name)(c); // call your kernel
}
};
@export(&Entry.run, .{ .name = name }); // name it for the host
}04 — the build
From Zig to a running shader, in three stages
So you've written one Zig file. How does it become a GPU shader? WebGPU in the browser doesn't run Zig, and it doesn't run the GPU's native machine code either — it runs WGSL, a high-level shading language. The journey has three steps:
- Zig → SPIR-V. Zig's own compiler can target SPIR-V, a portable binary format for GPU programs. This is just the normal Zig compiler with a GPU target.
- SPIR-V → WGSL. A custom tool in this project, spv2wgsl, reads that binary and writes out readable WGSL text. This is the part most engines outsource to a big C++ library; zimr does it in pure Zig.
- WGSL → GPU. The browser's WebGPU implementation compiles the WGSL the rest of the way and runs it.
The build wires the first two stages together. The Zig-to-SPIR-V step is an ordinary compiler invocation with a GPU target:
// Stage 1: Zig -> SPIR-V (same flags as the fragment-shader path).
const compile: *Run = b.addSystemCommand(&.{
b.graph.zig_exe, "build-obj",
"-target", "spirv32-vulkan",
"-mcpu", "vulkan_v1_2",
"-fno-llvm", "-fno-lld",
"-O", "ReleaseFast",
"-ofmt=spirv",
});
const raw_spv: LazyPath = compile.addPrefixedOutputFileArg("-femit-bin=", "compute.spv");Then the SPIR-V is handed to spv2wgsl, along with the workgroup size to stamp into the output:
// Stage 2: SPIR-V -> WGSL with the injected @workgroup_size.
w.addArg(b.fmt("--workgroup={d},{d},{d}", .{ workgroup[0], workgroup[1], workgroup[2] }));
if (entry) |e| { w.addArg(b.fmt("--entry={s}", .{e})); }
w.addFileArg(raw_spv);
return w.addOutputFileArg("compute.wgsl"); // the finished shader textThe resulting WGSL is baked into the app, where the host loads it onto the GPU. One source file, run through a pipeline, comes out the other end as a shader — with the matching CPU loop built from the very same source in a separate pass.
05 — the translator
spv2wgsl: turning GPU binary back into readable shaders
spv2wgsl is the most unusual piece. SPIR-V is a binary list of numbered instructions — closer to assembly than to source. WGSL is structured, high-level text with real ifs and fors. Going from the first to the second means rebuilding structure that the binary has flattened: recovering loops from jumps, naming anonymous values, mapping each SPIR-V operation to its WGSL equivalent. It's a small compiler in its own right, and in this project it's around ten thousand lines of Zig.
Most of it is mechanical translation. The interesting parts are where the two languages don't line up. Atomics are the headline example, and they show how the whole system absorbs a gap in the toolchain rather than giving up.
The atomics gap
The version of Zig used here can target SPIR-V, but its SPIR-V backend doesn't yet implement atomic operations — the very thing the parallel grid build depends on. You'd be stuck, except for a trick. Recall the k.atomicAdd helper from earlier: on the GPU it calls a deliberately useless function named zatomicAdd:
// MUST be `noinline` so the call survives into the SPIR-V for spv2wgsl
// to find. The body is never run on the GPU — it exists only to make
// Zig emit a real function call with the right name.
noinline fn zatomicAdd(arr: anytype, idx: u32, val: u32) u32 {
const prev: u32 = arr[idx];
arr[idx] = prev +% val;
return prev;
}Because it's marked noinline, Zig can't optimize it away — it leaves a real function call in the SPIR-V, named after the helper. spv2wgsl watches for that name. When it sees a call to something containing zatomicAdd, it does four things: it emits the genuine WGSL atomicAdd(&field[idx], val) instead of the call, it deletes the useless helper function, it marks the target array as atomic, and it changes that array's type to array<atomic<u32>> in the output. The placeholder is a hook; spv2wgsl rewrites it into the real builtin.
A one-character bug that froze the fluid
The translator also has to get small things exactly right, and one of them caused a memorable failure. In SPIR-V, a value can be given a name more than once. A data array that gets passed into a function picks up the function's parameter name as a second, later name. spv2wgsl was using "last name wins," so the array grid_counts ended up renamed to the generic parameter name arr. The host, which finds each buffer by its real name, then bound the wrong memory — and on the device the fluid froze into a dead blob. The fix is one line, and a comment now guards it: first name wins.
// Zig emits MULTIPLE names per id; a storage array passed as a function
// arg picks up the callee's PARAMETER name as a later duplicate. Keep the
// FIRST name so the host's `kbuf_` lookup still matches the buffer.
if (s.ids[target].wgsl_name.len == 0) s.ids[target].wgsl_name = name; The lesson worth taking from this section isn't the specific bug — it's that a translator between two languages is a sequence of decisions this precise, and each one can be the difference between a working simulation and a frozen screen.
06 — the host
The host: starting the work, and keeping data on the GPU
The host is the ordinary code that runs on the CPU and tells the GPU what to do: load the shaders, allocate the buffers, and launch kernels. The launch function, run, is where the CPU/GPU split appears one more time. On the CPU it's literally a loop; on the GPU it records a dispatch:
pub fn run(self: *Self, comptime name: []const u8) void {
switch (self.backend) {
.cpu => {
// the GPU's "thousands of threads" is just a for-loop here
var id: u32 = 0;
while (id < self.count) : (id += 1) {
@field(M, name)(.{ .id = id, .params = self.params });
}
},
.gpu => {
// ...find the compiled pipeline for `name`, then:
const wg: u32 = M.config.workgroup;
compute_pass.setPipeline(cp, pipeline);
compute_pass.setBindGroup(cp, 0, gp.bind_group);
// launch ceil(count / workgroup) groups of `workgroup` threads
compute_pass.dispatchWorkgroups(cp, .{ .x = (self.count + wg - 1) / wg });
},
}
}Notice dispatchWorkgroups launches count / workgroup groups — that's why the in-kernel guard matters: the division rounds up, so the last group has spare threads.
One submission per frame
Talking to the GPU has overhead. Every time the host "submits" work, the driver has paperwork to do. A naive simulation submits once per kernel — and with several kernels run several times per frame, that paperwork dominates. The fix is batching: open one recording, encode all the kernels and all the sub-steps into it, and submit once.
/// Open a single-pass batch: every following `run` encodes into one
/// compute pass until `endBatch` submits it.
pub fn beginBatch(self: *Self) void {
// ... write params once, open one encoder + one compute pass
gp.batch_pass = compute_pass.begin(gp.batch_enc);
compute_pass.setBindGroup(gp.batch_pass, 0, gp.bind_group);
}
pub fn endBatch(self: *Self) void {
compute_pass.end(gp.batch_pass);
const cmd = wgpu.finishCommandEncoder(gp.batch_enc);
wgpu.queueSubmit(gp.queue, cmd); // ONE submit for the entire frame
}With batching on, the kernels just record themselves into the open pass and nothing is submitted until the frame is fully assembled.
The data never leaves the GPU
The most important performance fact about this fluid is something the host doesn't do. Once the particles are on the GPU, they stay there. The compute kernels write positions into a GPU buffer, and the renderer reads positions from that same buffer to draw the discs — no copy back to the CPU, no round trip. The simulation and the drawing share memory.
The only time data is copied back to the CPU is for the on-screen diagnostics (the "cells / average" readout), and that copy is switched off unless you open the diagnostics panel. The function that does it, readLatest, also carries a small scar:
// Variable-count fields (positions, …) are sized to the particle count;
// fixed-size fields (grid_counts, cell_start) must slice to THEIR OWN
// length, never the particle count — otherwise grid_counts overruns into
// the cell_start array that follows it. Clamp to the lesser of the two.
const arr_len: usize = @field(M.g.B, @tagName(field)).len;
const n: usize = @min(self.count, arr_len);Before this fix, the diagnostics asked for 20,000 grid cells when only ~1,300 existed, read straight past the end into the next array, and reported impossible numbers (a "fullest cell" of 19,999). The simulation itself was always fine — only the readout was lying. It's a clean illustration of how a bug in a measurement can look exactly like a bug in the thing being measured.
07 — the example
The fluid, kernel by kernel
Now the payoff: the actual simulation. The method is Clavet's double-density relaxation, a particle fluid that's stable and cheap enough for real time. Each frame runs a few sub-steps, and each sub-step is a short sequence of kernels. First, the data they all share:
pub const num_particles: u32 = 20000;
pub const interact_radius: f32 = 22.0; // a particle only feels neighbours within this
pub const Buffers = extern struct {
pos: [num_particles]zm.Vec2, // position (pos is FIRST: the renderer reads it)
prev: [num_particles]zm.Vec2, // last position (used to recover velocity)
vel: [num_particles]zm.Vec2, // velocity
delta: [num_particles]zm.Vec2, // pending position correction
density: [num_particles]zm.Vec2, // (density, near-density) — also drives colour
grid_counts: [grid_cells]u32, // particles per cell (atomic); reused as the scatter cursor
cell_start: [grid_cells + 1]u32, // prefix sum: cell c owns sorted slots [cell_start[c] .. cell_start[c+1])
pos2: [num_particles]zm.Vec2, // scatter scratch — the cell-sorted copy of pos/vel/prev
vel2: [num_particles]zm.Vec2,
prev2: [num_particles]zm.Vec2,
};The neighbour grid: a counting sort, so the reads are fast
Checking all 20,000 particles against each other would be 400 million comparisons. Instead the domain is chopped into a grid of square cells, each the size of the interaction radius. A particle only needs to look in its own cell and the eight around it — a 3×3 neighbourhood — because anything farther away is out of reach.
The interesting question is how a cell's particles are stored. The obvious way is a list of indices per cell; the neighbour walk then reads pos[index] for each one — and those indices point all over memory, so every read is a random jump. On a GPU that is the slowest thing you can do (see the performance section). The fast way is to physically reorder the particles into cell order every sub-step, so that a cell's particles occupy one contiguous run of memory. That reordering is a classic counting sort, and it's five small kernels.
clearGrid zeroes every cell's counter — one thread per cell:
pub fn clearGrid(c: k.Ctx(@This())) void {
const cell: u32 = c.id;
if (cell >= c.params.n_cols * c.params.n_rows) return;
k.atomicStore(b_grid_counts, cell, 0);
}countGrid tallies how many particles fall in each cell — one thread per particle, an atomic add so thousands can count into the same cell at once:
pub fn countGrid(c: k.Ctx(@This())) void {
const i: u32 = c.id;
if (i >= c.params.count) return;
const cell: u32 = cellOf(b_pos[i], c.params.n_cols, c.params.n_rows, c.params.h);
_ = k.atomicAdd(b_grid_counts, cell, 1); // result unused — we only want the totals
}prefixSum turns those per-cell counts into start offsets. A running total: cell 0 starts at 0, cell 1 starts after cell 0's particles, and so on. After it, cell c owns the contiguous slot range [cell_start[c], cell_start[c+1]). It's a single serial scan (~1,300 cells — nothing next to 20,000 particles), so it runs on one thread, and it zeroes grid_counts as it goes so the next kernel can reuse it as a write cursor:
pub fn prefixSum(c: k.Ctx(@This())) void {
if (c.id != 0) return; // a single thread does the whole scan
const cells: u32 = c.params.n_cols * c.params.n_rows;
var acc: u32 = 0;
var cell: u32 = 0;
while (cell < cells) : (cell += 1) {
const cnt: u32 = k.atomicLoad(b_grid_counts, cell);
b_cell_start[cell] = acc; // where this cell's run begins
acc += cnt;
k.atomicStore(b_grid_counts, cell, 0); // reset → reused as the scatter cursor
}
b_cell_start[cells] = acc; // == num_particles
}scatter does the actual move: each particle claims the next free slot in its cell (atomic add on the reused counter) and writes its position, velocity, and previous-position into the sorted scratch at that slot. copyback then copies the scratch back over the originals, so pos/vel/prev become the cell-sorted order — the canonical layout for the rest of the sub-step and the next one. Because particles are interchangeable, there's no need to ever un-sort:
pub fn scatter(c: k.Ctx(@This())) void {
const i: u32 = c.id;
if (i >= c.params.count) return;
const cell: u32 = cellOf(b_pos[i], c.params.n_cols, c.params.n_rows, c.params.h);
const slot: u32 = k.atomicAdd(b_grid_counts, cell, 1); // unique offset within the cell
const dest: u32 = b_cell_start[cell] + slot; // its place in the sorted array
b_pos2[dest] = b_pos[i];
b_vel2[dest] = b_vel[i];
b_prev2[dest] = b_prev[i];
}There's a quiet bonus: with exact ranges there's no fixed per-cell capacity to overflow, so dense clumps never silently drop neighbours the way a capped slot list could.
Density: how crowded am I?
density walks the 3×3 neighbourhood and sums up how close the neighbours are. Two numbers come out: a normal density and a sharper "near-density" that spikes when particles are almost touching. The near-density is what stops the fluid from collapsing into itself. Because of the sort, each neighbour cell is just a contiguous slice — the loop reads consecutive particles from consecutive memory:
// for each of the 9 neighbour cells: walk its contiguous run
var kk: u32 = b_cell_start[cell];
while (kk < b_cell_start[cell + 1]) : (kk += 1) { // consecutive kk → consecutive memory
if (kk == i) continue;
const sep: zm.Vec2 = b_pos[kk] - my_pos;
const d2: f32 = zm.dot(sep, sep);
if (d2 < c.params.h * c.params.h) { // squared compare: skip the sqrt if out of range
const dist: f32 = @sqrt(d2);
const omq: f32 = 1.0 - dist / c.params.h;
rho += omq * omq; // density
rho_near += omq * omq * omq; // near-density (sharper)
}
}Force: relieve the pressure
force turns those densities into movement. A crowded particle pushes its neighbours away; the push is stronger the closer they are. Rather than write the pushes straight to velocity, the kernel accumulates a position correction, which the final kernel turns back into velocity. It reads the same contiguous neighbour runs the density pass did — now also pulling each neighbour's density and velocity from consecutive memory:
const omq: f32 = 1.0 - dist / c.params.h;
const jd: zm.Vec2 = b_density[kk]; // neighbour density — contiguous read
const disp: f32 = 0.5 * c.params.dt * c.params.dt *
((my_press + j_press) * omq + (my_near + j_near) * omq * omq);
corr = corr - dir * zm.splat2(disp); // push away from the crowdViscosity: a separate, gentler pass
viscosity is what makes the fluid look like water rather than a gas — a drag between particles moving toward each other. It used to be folded into the force loop, but it now runs as its own pass before the prediction step, matching the order in Clavet's paper, which is noticeably more stable. The model is a small collision impulse: it acts only on a pair that's approaching (closing speed u > 0), it's quadratic in that speed, and it's clamped so it can at most bring the pair's approach to rest — never a bounce, so it can only ever calm the fluid, never inject energy:
const u: f32 = zm.dot(my_vel - b_vel[kk], n); // closing speed along the line of centres
if (u > 0.0) { // only if approaching
const w: f32 = if (w_lin < 1.0) w_lin else 1.0; // 1 inside h/2, fading to 0 at h
const imp_raw: f32 = c.params.visc_beta * w * u * u; // quadratic in the closing speed
const imp: f32 = if (imp_raw < u) imp_raw else u; // clamp: never overshoot to a bounce
my_vel = my_vel - n * zm.splat2(0.5 * imp); // half to each side
}Gravity, predict, finalize: move, and hit the walls gracefully
Three per-particle kernels bracket the neighbour work. gravityMouse applies gravity and any mouse push to the velocity. predict remembers each particle's current position (into prev) and advances it by its velocity. applyAndFinalize (run last) adds the accumulated correction, keeps particles inside the box, and recovers each velocity from how far it actually moved.
The wall handling here fixed a bug that resisted fixes for a long time: particles exploding in the corners. The cause was a hard clamp that pinned a cornered particle to the exact corner point on both axes at once. Several particles landing on the identical point made the near-density term blow up — and they shot apart. The fix is a soft push that ramps up near a wall (so density stays smooth), plus a tiny per-particle nudge along the wall so two particles never share one point:
const band: f32 = c.params.h; // start pushing one radius from the wall
const push: f32 = 0.25;
if (p[0] < band) { p[0] += (band - p[0]) * push; }
else if (p[0] > c.params.dom_w - band) { p[0] -= (p[0] - (c.params.dom_w - band)) * push; }
// ... same for the vertical walls ...
// hard safety clamp, with a sub-particle jitter ALONG the wall so two
// cornered particles never land on the exact same point (which spiked density)
const jitter: f32 = @as(f32, @floatFromInt(i % 7)) * 0.3; // 0..1.8 pxThe whole sub-step, on one line
Put together, each sub-step the host records is eleven dispatches — three per-particle steps, the five-kernel sort, then the two neighbour walks and the finalize:
s.pipe.run("gravityMouse"); // gravity + mouse → velocity
s.pipe.run("viscosity"); // collision-impulse drag (on last step's sorted grid)
s.pipe.run("predict"); // remember prev, advance by velocity
s.pipe.run("clearGrid"); // ─┐
s.pipe.run("countGrid"); // │ counting sort: reorder pos/vel/prev
s.pipe.run("prefixSum"); // │ into cell order so the neighbour
s.pipe.run("scatter"); // │ walks below read contiguous memory
s.pipe.run("copyback"); // ─┘
s.pipe.run("density"); // measure crowding (neighbour walk #1)
s.pipe.run("force"); // push apart (neighbour walk #2)
s.pipe.run("applyAndFinalize"); // apply, clamp to walls, recover velocityThat sequence, run three times per frame for 20,000 particles, batched into a single submit, with the results drawn straight from GPU memory — that's the whole fluid. (Viscosity runs before the sort, so it uses last sub-step's ordering — a damping term doesn't mind a slightly stale grid; density and force run after, on the fresh one.)
08 — performance
Making it fast: what every optimization actually does
A simulation can be correct and still too slow to use. The number that matters here is milliseconds per sub-step — how long one of those kernel sequences takes on the phone. Lower is better; to hit a smooth 60 frames per second with three sub-steps per frame, the budget is tight. Here's the path the project actually walked, from where it started to where it is now:
At three sub-steps a frame, 3.8 ms each comes to about 11 milliseconds of compute — so the simulation clears the 60-frames-per-second bar with real headroom, and on the test phone (an Adreno GPU) it runs at close to 90. That's the whole point of measuring in milliseconds per sub-step rather than frames per second: the frame rate hides behind the display's refresh cap, but the sub-step time is the honest throughput number, and it's what dropped.
The last step is the big one — it nearly tripled the throughput in a single change. It's also the most important idea in this whole document, so it's worth being precise about why it works.
The headline: coalesced memory, and why the sort pays for itself
On a GPU, threads run in lockstep groups (a "warp" — 32 or 64 of them). When all the threads in a group read memory at the same instant, the hardware can satisfy them efficiently only if the addresses they want are next to each other — then one wide memory transaction serves the whole group. That's called a coalesced read. If the addresses are scattered, the hardware must issue a separate transaction per thread, and the group stalls while they trickle in one by one. Same number of bytes, many times the latency.
Now look at the neighbour walk. The slow grid stored a list of particle indices per cell, so reading a neighbour meant b_pos[grid_data[slot]] — an index pulled from one place, then a position fetched from wherever that particle happened to live in memory. Adjacent threads, looking at the same cell, chased indices pointing all over the buffer. Every neighbour read was a scattered, uncoalesced jump. For a fluid, the neighbour walks are nearly the entire cost, so this one access pattern dominated the frame.
The counting sort fixes it at the source: after the sort, a cell's particles sit in one contiguous run of the pos array. The neighbour loop just walks b_pos[kk] for consecutive kk — consecutive addresses — so the reads coalesce, and the same is true for the density and velocity the force pass reads. The fluid does the identical math; it just stops fighting the memory system. The sort itself costs four extra little passes (count, the one-thread scan, scatter, copy-back), but every one of them is cheap and coalesced, and they buy back far more than they cost in the two neighbour walks that follow.
Three earlier ideas are still in the code and still pull their weight:
- Atomics for parallel binning.
countGridandscatteruse atomic adds so all 20,000 particles can tally into — and claim slots in — shared cells at once, with the hardware guaranteeing no two collide. Atomics are the price of admission for "many threads writing to shared structures safely." - Workgroup size 256. The GPU runs threads in groups; too small and the chip sits half-idle waiting on memory with nothing else to do. Launching in groups of 256 (up from 64) gives the scheduler enough threads that while some wait for data, others compute — exactly what hides the memory latency the sort doesn't eliminate.
- Skip the square root. To tell if a neighbour is in range you compare squared distances, avoiding a square root for the many candidates that turn out too far away — the 3×3 box reaches past the circular radius, so its corners are always out of range and now cost nothing.
09 — the road ahead
The turn we didn't take, and what's actually left
For a long time the plan in this very section was shared-memory tiling: assign one thread-group to one cell, have it copy the cell's neighbourhood onto the GPU's small, fast on-chip scratchpad once, then let every particle in the cell read neighbours from the scratchpad instead of making repeated trips to main memory. On paper it's the textbook fix for redundant neighbour reads, and it's the kind of thing that closes the gap to a hand-tuned reference. So it got built — and on this fluid it was about 19% slower.
Why tiling lost
The killer was occupancy. Tiling wants one workgroup per cell, and a workgroup here is 256 threads — but a typical cell holds only ~17–30 particles. So fewer than one lane in eight had a particle to work on; the rest sat idle holding the group open. The scratchpad did save memory traffic, but the machine was so under-filled that it lost far more than it saved. A clean illustration that an optimization can be locally correct — fewer reads! — and globally a loss.
Why the sort won instead
The counting sort chases the same prize — make the neighbour reads cheap — but from the opposite direction. Instead of caching scattered data in fast memory, it rearranges the data so it isn't scattered in the first place, and it keeps the simple, fully-occupied "one thread per particle" launch that fills the machine. Same goal (less effective memory traffic), no occupancy penalty, and a much simpler kernel: no scratchpad, no barrier. It's the better lever, and it's the one that actually moved the number — 10 ms to 3.8.
What's actually left
With the big lever pulled, the rest is tuning, not restructuring: choosing a particle count and cell size that hold a steady 60 frames per second with headroom, and trimming the cost of drawing 20,000 overlapping discs. One structural option remains in reserve — the prefix sum is a single-thread serial scan today, which is fine at ~1,300 cells but would become a parallel scan if the grid ever grew large enough for it to show up in a profile. It doesn't yet. Smaller knobs, turned now that the fluid is fast.