WebGPU notes   /   37

WebGPU N-Body Simulation in Rust: 12,288 GPU Particles

Stage 256 source positions at a time in workgroup memory, evaluate all 150,994,944 particle pairs, integrate 12,288 bodies in place, and expose the attraction and billboard parameters through a live egui panel.

By Pooya Eimandar·Rust / wgpu / WGSL

Particles

12,288

Pair terms per frame

150,994,944

Shared tile per workgroup

4 KiB

Tile an all-pairs N-body calculation in workgroup memory

The previous compute culling and LOD example uses compute work to prepare only the mesh draws needed for a camera view. Compute N-body applies the same programmable stage to simulation instead. Every body contributes attraction to every target, and one storage buffer remains in place. A calculate pass reads positions and writes velocities; a second compute pass integrates those new velocities into positions. The render pass then reads the same buffer as vertex-stage storage.

A direct all-pairs shader would issue one global position load for each of 12,288² interactions. This shader instead has each 256-thread workgroup load one 256-position tile into 4 KiB of shared memory. All 256 target threads reuse that tile before advancing. Tiling reduces source-position storage loads by a factor of 256 while preserving the full quadratic arithmetic.

The introducing Compute N-body commit added the Rust source, compute WGSL, render WGSL, screenshots, build registration, README entry, and gallery registration. It reused the existing particle textures. Current code later received the modern WebAssembly entry, direct WebGPU projection, seed-shape and finite-distance fixes, and collapsible settings window. The shared loader's later batch update supplies the current concurrent browser fetch. Rendering runs through sib::render.

Load two uncompressed KTX 1 textures

Before entering sib::render, the shared asset loader requests particle01_rgba.ktx and particle_gradient_rgba.ktx as one batch. Native gives each request a thread. The current WebAssembly path sends both URLs to one temporary Worker and fetches them concurrently with Promise.all.

The example's local decoder accepts little-endian KTX 1 data with uncompressed unsigned-byte RGBA8, one face, no array, and no volume depth. It converts a zero image height to one, copies every encoded mip into a CPU RGBA vector, and uploads Rgba8Unorm. Both inputs have one mip, so the configured linear mip filter has no lower level to select.

Scroll sideways to see all table columns.

The two runtime files total 17,608 stored bytes and upload 17,408 texel bytes. Both are linear RGBA8 rather than sRGB, and both alpha channels are fully opaque. They are exact byte matches for the same particle sprite and gradient in Sascha Willems' Vulkan-Assets repository. Its README gives no general asset or per-file license, so this article does not assign definitive authorship or licensing to the textures.

Seed six clusters with heavy and ordinary bodies

A deterministic linear congruential generator starts from 0x4d2f_6a31. Rust creates 2,048 particles around each of six seed locations: ±5 on X, ±5 on Z, Y = 4, and Y = −8. The first record in every cluster begins at 1.5 times its seed location with mass 90,000. The remaining 2,047 first choose an exact random-direction 0.75 offset, then multiply absolute world-space Y by 2 - 0.75² = 1.4375. Their final positions are therefore not on a 0.75-radius shell, and the two Y-axis clouds shift relative to their nominal seeds. They receive masses from 37.5 to 75 and tangential velocities plus jitter; rotation direction alternates between clusters.

These six heavy bodies are not pinned attractors. They occupy the same buffer, contribute to every target, receive forces, and integrate like the other 12,282 records. The cluster index divided by six initializes each body's gradient phase in vel.w.

Scroll sideways to see all table columns.

The fixed total excludes the 4 KiB shared array allocated per resident calculate workgroup, bind-group and sampler descriptors, presentation surface, egui resources, embedded font, staging behavior, and driver alignment. “Reset particles” regenerates the deterministic CPU vector and writes all 393,216 bytes back into the storage buffer without reallocating it.

Calculate 150,994,944 pair terms with 48 tiled workgroups

The calculate pass dispatches 48 workgroups of 256 threads, exactly covering all 12,288 targets. Every workgroup walks 48 source tiles. Each lane loads one vec4 containing XYZ and mass, all lanes synchronize, and each target accumulates 256 attraction terms from shared memory. A second barrier prevents the next tile from overwriting positions while another lane still reads them.

for (var tile_start = 0u;
     tile_start < particle_count;
     tile_start = tile_start + 256u) {
    let source_index = tile_start + local_id.x;
    if (source_index < particle_count) {
        shared_positions[local_id.x] = particles[source_index].pos;
    } else {
        shared_positions[local_id.x] = vec4<f32>(0.0);
    }
    workgroupBarrier();

    for (var i = 0u; i < 256u; i = i + 1u) {
        let other = shared_positions[i];
        let delta = other.xyz - position.xyz;
        let dist_sq = max(dot(delta, delta) + soften, 1.0e-6);
        acceleration += gravity * delta * other.w / pow(dist_sq, power);
    }
    workgroupBarrier();
}

Scroll sideways to see all table columns.

The 9 MiB figure counts tiled source-position loads only. Adding initial target position and velocity reads, calculate velocity stores, integrate position and velocity reads, and integrate position stores brings the two compute shaders to 10,616,832 logical storage bytes, or 10.125 MiB, per frame. Of that, 393,216 bytes are shader storage writes. These are source-level counts, not measured memory traffic or execution time. Caches, transactions, compiler transformations, occupancy, and the cost of 151 million dynamic pow evaluations determine physical performance. Tiling attacks global reads; it does not reduce O(N²) force arithmetic.

The force is an adjustable visual law: gravity * delta * mass / (distance_squared + soften)^power. Softening is added directly to squared distance rather than being squared itself. With default power 0.75 the law is not Newton's inverse-square acceleration, which would require an exponent of 1.5 when multiplying by the unnormalized delta; the UI stops at 1.4. Below power 0.5, including part of the exposed range, unsoftened force magnitude increases with distance. The 1e−6 floor added in the later fix prevents a zero-over-zero self term from becoming NaN even if softening is set outside the UI range.

Integrate updated velocity in a second compute pass

The calculate pass writes only vel, while all workgroups continue to read pos. That field separation makes one in-place buffer safe without ping-ponging. A distinct integrate pass provides global ordering, then adds the updated velocity to position using the same time step. This is a split, semi-implicit Euler update.

Current WGSL adds the full velocity vec4 to the full position vec4. XYZ behaves as expected, but vel.w is the animated gradient phase and pos.w is mass. The phase therefore also increases mass by delta_t * vel.w every active frame. At default timing, an average phase near 0.5 adds roughly 0.025 mass units per real second. Because mass controls both attraction strength and billboard size, the simulation slowly changes those properties even though Rust initialized them as fixed values. Resetting particles restores the original masses.

The uniform time step is zero while paused. Otherwise it is clamp(frame_delta * 0.05 * time_scale, 0, 1/120). At the default time scale, simulation time advances at 5% of wall time, and the gradient phase advances about 0.005 per real second for a roughly 200-second cycle. Pause and time scale zero freeze numerical state, but both compute passes still dispatch and the calculate pass still evaluates every pair.

Render 12,288 projected textured billboards

The graphics pipeline has no vertex or index buffer. One instanced call emits six procedural corners per particle: 73,728 vertex invocations and 24,576 triangles. The vertex shader reads position, mass, and phase from storage. A fixed 60° perspective camera sits 14 units from the origin at 75° yaw and 26° pitch, with near and far planes of 0.1 and 512.

Projected mass and depth produce a physical point size, multiplied by the Particle size control and clamped from 1 to 128 surface pixels. The six 90,000-mass bodies reach the upper clamp. Each quad remains screen-aligned. There is no depth attachment, depth test, culling, or multisampling, so every covered fragment blends regardless of which system is nearer.

The fragment shader multiplies the radial sprite by the phase-selected gradient and brightness. RGB uses additive one-plus-one blending, so alpha would not gate RGB even if the sprite contained variable alpha. Alpha uses source one and destination one-minus-source-alpha. Because the KTX sprite alpha is opaque everywhere and the surface clears with alpha one, final surface alpha remains one. Dark sprite corners still execute and add near-black RGB because the shader has no discard. With every body on screen at the 128×128 size clamp, the loose upper bound is 201,326,592 fragment candidates before clipping and subpixel rasterization.

Tune seven controls without rebuilding GPU resources

The live egui panel is fixed-width, non-resizable, and collapsible. Changing a control rewrites the same 176-byte uniform without rebuilding pipelines, bindings, textures, or buffers. “Reset params” restores every default, including unpausing. “Reset particles” preserves controls and replaces particle state with the deterministic initial vector.

Scroll sideways to see all table columns.

egui owns input; there is no camera, pointer-force, or keyboard control outside the panel. Its displayed frame time and FPS average CPU or event-loop delivery rather than GPU timestamps. The screenshot shows time scale 2.20, power 0.80, soften 0.070, and brightness 2.45, not the defaults above.

Record two compute passes and two render passes

The calculate and integrate entry points run in separate compute passes. A particle pass clears the presentation surface and draws every body, then the egui pass loads that color target and renders the settings panel. All graphics work is single-sample.

Scroll sideways to see all table columns.

A controls-stable frame performs one 176-byte example uniform queue write before any egui uploads. Projection and the fixed model-view matrix are recalculated and included every frame even though only aspect ratio can change the projection. A control change causes another 176-byte uniform write during GUI rendering; because it is queued before the main encoder submission, the later values affect that submitted frame. Reset particles adds one 393,216-byte buffer write before submission, so the compute passes can immediately advance the reset state once unless paused. The renderer records the same four passes in all three cases.

Run and extend the example

From a local checkout with Rust installed, run the native WebGPU Compute N-body example:

cargo run --example computenbody

For WebAssembly, install the wasm32-unknown-unknown target and the wasm-bindgen CLI version matching Cargo.lock, then build and serve:

scripts/build-wasm.sh --release computenbody
cargo run --bin serve

Open http://127.0.0.1:8080/computenbody/ in a browser with WebGPU support. Both targets load the same two KTX files before calling sib::render, then run the same tiled calculate, integrate, billboard, and GUI passes.

The algorithm remains O(N²). Tiling saves global reads but does not reduce interaction count, and every pair uses a general pow. There is no Barnes–Hut tree, fast multipole method, neighbor cutoff, adaptive body count, subgroup optimization, GPU timing, or capability-based workload selection. Doubling particle count would approximately quadruple pair arithmetic. The fixed 256 appears in Rust dispatch math, both WGSL workgroup attributes, the shared-array length, and the tile step; those values must stay synchronized.

The force exponent is artistic rather than Newtonian, and single-precision semi-implicit Euler integration has no adaptive step or energy correction. Force terms accumulate in a fixed tile order without compensated summation. Bodies do not collide, merge, conserve a fixed system energy, or leave trails in a history buffer. Close approaches depend strongly on gravity, exponent, softening, and the capped time step. Pausing does not skip the expensive calculate pass.

The full-vec4 integrate operation couples gradient phase into mass and projected size. That drift is small for the six 90,000-mass bodies but meaningful over time for ordinary 37.5–75 masses. Updating only pos.xyz, keeping mass in a read-only field, or separating render phase would preserve the intended quantities.

The fixed camera cannot orbit, pan, or zoom. Additive billboards have no depth test, sorting, culling, indirect draw, soft-particle intersection, HDR target, exposure, tone mapping, or MSAA. Their opaque sprite alpha means black corners still execute fragment work. One-mip linear textures and a non-sRGB gradient further favor a compact demo over a calibrated rendering pipeline.

Useful changes to try:

• Update only position XYZ, keep mass fixed, and compare trajectories before and after the correction.
• Replace all-pairs evaluation with Barnes–Hut or a hierarchical approximation, then chart error against frame time.
• Use an approximate inverse-power path for common exponents and profile it against dynamic pow.
• Add camera controls, depth-aware billboards, HDR composition, GPU timestamps, and a runtime particle-count selector.