Skip to content
PNGine

Fluid simulation

Download PNG

A small stable fluid. Two fields live on a 96×96 grid, a velocity (vec2f per cell) and a dye density (f32 per cell); each frame a compute shader advects both by looking backwards along the velocity and sampling with bilinear interpolation, then stirs in two rotating sources. The display colours the flow’s direction as a hue and its speed as brightness, with the dye added on top. The poster is the first frame; press Play to see it swirl.

examples/samples/19_fluid_simulation.sjon
; Fluid simulation. A 96×96 semi-Lagrangian fluid: each frame a compute pass
; traces velocity and density backward along the flow with bilinear sampling,
; adds two rotating force/dye sources and damps the boundary; a fullscreen
; triangle then colours velocity direction (HSV hue, speed as value) over
; density brightness. Two independent `:pool 2` ping-pong fields (velocity
; vec2f, density f32) give the step bind group five entries (uniforms + vel
; in/out + den in/out with `:ping-pong`); `(dispatch :workgroups [6 6 1])` over a 16×16
; workgroup. A compute `(init …)` zeroes the velocity field once.
;
; GRID_SIZE sizes the buffers only; each shader restates the extent as a local
; `SIZE` constant. Keep them in step.


(define :name GRID_SIZE :value 96)


(buffer :name velocityBuffer :size (* GRID_SIZE GRID_SIZE 2 4) :usage [storage] :pool 2)


(buffer :name densityBuffer :size (* GRID_SIZE GRID_SIZE 4) :usage [storage] :pool 2)


(buffer :name uniforms :size 16 :usage [uniform copy-dst])


(queue :name writeUniforms
  (write-buffer :buffer uniforms :offset 0 :data pngine-inputs))


; Zero the velocity field once.
(init :name initFluid :buffer velocityBuffer :module initShader :workgroups [(ceil (/ (* GRID_SIZE GRID_SIZE) 64))])


(shader-module :name initShader :code """
struct Velocity { data: array<vec2f> }
struct Density { data: array<f32> }


@binding(0) @group(0) var<storage, read_write> velocity: Velocity;


const SIZE: u32 = 96u;


@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) id: vec3u) {
  let idx = id.x;
  if (idx >= SIZE * SIZE) { return; }
  velocity.data[idx] = vec2f(0.0);
}
""")


(shader-module :name stepShader :code """
struct Uniforms {
  time: f32,
  width: f32,
  height: f32,
  aspect: f32,
}
@group(0) @binding(0) var<uniform> u: Uniforms;


struct Velocity { data: array<vec2f> }
struct Density { data: array<f32> }


@group(0) @binding(1) var<storage, read> velIn: Velocity;
@group(0) @binding(2) var<storage, read_write> velOut: Velocity;
@group(0) @binding(3) var<storage, read> denIn: Density;
@group(0) @binding(4) var<storage, read_write> denOut: Density;


const SIZE: u32 = 96u;
const SIZEF: f32 = 96.0;
const DT: f32 = 0.5;
const DECAY: f32 = 0.995;


fn sampleVelBilinear(x: f32, y: f32) -> vec2f {
  let x0 = clamp(i32(floor(x)), 0, i32(SIZE) - 1);
  let y0 = clamp(i32(floor(y)), 0, i32(SIZE) - 1);
  let x1 = min(x0 + 1, i32(SIZE) - 1);
  let y1 = min(y0 + 1, i32(SIZE) - 1);


  let fx = fract(x);
  let fy = fract(y);


  let v00 = velIn.data[u32(y0) * SIZE + u32(x0)];
  let v10 = velIn.data[u32(y0) * SIZE + u32(x1)];
  let v01 = velIn.data[u32(y1) * SIZE + u32(x0)];
  let v11 = velIn.data[u32(y1) * SIZE + u32(x1)];


  return mix(mix(v00, v10, fx), mix(v01, v11, fx), fy);
}


fn sampleDenBilinear(x: f32, y: f32) -> f32 {
  let x0 = clamp(i32(floor(x)), 0, i32(SIZE) - 1);
  let y0 = clamp(i32(floor(y)), 0, i32(SIZE) - 1);
  let x1 = min(x0 + 1, i32(SIZE) - 1);
  let y1 = min(y0 + 1, i32(SIZE) - 1);


  let fx = fract(x);
  let fy = fract(y);


  let d00 = denIn.data[u32(y0) * SIZE + u32(x0)];
  let d10 = denIn.data[u32(y0) * SIZE + u32(x1)];
  let d01 = denIn.data[u32(y1) * SIZE + u32(x0)];
  let d11 = denIn.data[u32(y1) * SIZE + u32(x1)];


  return mix(mix(d00, d10, fx), mix(d01, d11, fx), fy);
}


@compute @workgroup_size(16, 16)
fn main(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= SIZE || id.y >= SIZE) { return; }


  let x = f32(id.x);
  let y = f32(id.y);
  let idx = id.y * SIZE + id.x;


  let vel = velIn.data[idx];


  // Semi-Lagrangian advection: trace backward
  let prevX = x - vel.x * DT;
  let prevY = y - vel.y * DT;


  var newVel = sampleVelBilinear(prevX, prevY);
  var newDen = sampleDenBilinear(prevX, prevY);


  // Add rotating force sources
  let t = u.time;


  // Source 1 - red
  let cx1 = SIZEF * 0.5 + sin(t * 0.7) * SIZEF * 0.3;
  let cy1 = SIZEF * 0.5 + cos(t * 0.5) * SIZEF * 0.3;
  let d1 = length(vec2f(x - cx1, y - cy1));
  if (d1 < 6.0) {
    let force = vec2f(cos(t * 2.0), sin(t * 2.0)) * 8.0;
    newVel += force * (1.0 - d1 / 6.0);
    newDen += 2.0 * (1.0 - d1 / 6.0);
  }


  // Source 2 - blue
  let cx2 = SIZEF * 0.5 + sin(t * 0.9 + 2.0) * SIZEF * 0.25;
  let cy2 = SIZEF * 0.5 + cos(t * 0.6 + 1.0) * SIZEF * 0.25;
  let d2 = length(vec2f(x - cx2, y - cy2));
  if (d2 < 5.0) {
    let force = vec2f(sin(t * 1.5), -cos(t * 1.5)) * 6.0;
    newVel += force * (1.0 - d2 / 5.0);
    newDen += 1.5 * (1.0 - d2 / 5.0);
  }


  // Boundary conditions
  if (id.x < 2u || id.x >= SIZE - 2u || id.y < 2u || id.y >= SIZE - 2u) {
    newVel *= 0.5;
  }


  velOut.data[idx] = newVel * DECAY;
  denOut.data[idx] = clamp(newDen * DECAY, 0.0, 3.0);
}
""")


(compute-pipeline :name stepPipeline :layout auto (compute :module stepShader :entry main))


(bind-group :name stepBindGroup :layout stepPipeline :group 0 :pool 2
  (entry :binding 0 :buffer uniforms)
  (entry :binding 1 :buffer velocityBuffer :ping-pong 0)
  (entry :binding 2 :buffer velocityBuffer :ping-pong 1)
  (entry :binding 3 :buffer densityBuffer :ping-pong 0)
  (entry :binding 4 :buffer densityBuffer :ping-pong 1))


(compute-pass :name stepPass
  :pipeline stepPipeline
  :bind-groups [stepBindGroup]
  :bind-groups-pool-offsets [0]
  (dispatch :workgroups [6 6 1]))


(shader-module :name renderShader :code """
struct Uniforms {
  time: f32,
  width: f32,
  height: f32,
  aspect: f32,
}
@group(0) @binding(0) var<uniform> u: Uniforms;


struct Density { data: array<f32> }
struct Velocity { data: array<vec2f> }


@group(0) @binding(1) var<storage, read> density: Density;
@group(0) @binding(2) var<storage, read> velocity: Velocity;


const SIZE: f32 = 96.0;


@vertex
fn vs(@builtin(vertex_index) i: u32) -> @builtin(position) vec4f {
  let x = f32(i & 1u) * 4.0 - 1.0;
  let y = f32((i >> 1u) & 1u) * 4.0 - 1.0;
  return vec4f(x, y, 0.0, 1.0);
}


fn hsv2rgb(h: f32, s: f32, v: f32) -> vec3f {
  let c = v * s;
  let x = c * (1.0 - abs(fract(h * 6.0) * 2.0 - 1.0));
  let m = v - c;


  var rgb = vec3f(0.0);
  let hi = u32(h * 6.0) % 6u;
  if (hi == 0u) { rgb = vec3f(c, x, 0.0); }
  else if (hi == 1u) { rgb = vec3f(x, c, 0.0); }
  else if (hi == 2u) { rgb = vec3f(0.0, c, x); }
  else if (hi == 3u) { rgb = vec3f(0.0, x, c); }
  else if (hi == 4u) { rgb = vec3f(x, 0.0, c); }
  else { rgb = vec3f(c, 0.0, x); }


  return rgb + vec3f(m);
}


@fragment
fn fs(@builtin(position) pos: vec4f) -> @location(0) vec4f {
  let cellX = u32(pos.x * SIZE / u.width);
  let cellY = u32(pos.y * SIZE / u.height);


  if (cellX >= u32(SIZE) || cellY >= u32(SIZE)) {
    return vec4f(0.0, 0.0, 0.0, 1.0);
  }


  let idx = cellY * u32(SIZE) + cellX;


  let den = density.data[idx];
  let vel = velocity.data[idx];


  // Color from velocity direction
  let speed = length(vel);
  let hue = atan2(vel.y, vel.x) / 6.283 + 0.5;


  let velocityColor = hsv2rgb(hue, 0.8, min(speed * 0.15, 1.0));


  // Combine velocity color with density brightness
  let densityBrightness = clamp(den * 0.4, 0.0, 1.0);


  let color = velocityColor * 0.6 + vec3f(densityBrightness);


  return vec4f(color, 1.0);
}
""")


(render-pipeline :name renderPipeline
  :layout auto
  (vertex :module renderShader :entry vs)
  (fragment :module renderShader :entry fs
    (target :format preferred-canvas-format)))


(bind-group :name renderBindGroup :layout renderPipeline :group 0 :pool 2
  (entry :binding 0 :buffer uniforms)
  (entry :binding 1 :buffer densityBuffer :ping-pong 1)
  (entry :binding 2 :buffer velocityBuffer :ping-pong 1))


(render-pass :name drawPass
  (color-attachment :view context-current-texture :clear-value [0 0 0 1] :load-op clear :store-op store)
  :pipeline renderPipeline
  :bind-groups [renderBindGroup]
  :bind-groups-pool-offsets [0]
  (draw :vertex-count 3))


(frame :name main
  :init [initFluid]
  :perform [writeUniforms stepPass drawPass])

examples/samples/19_fluid_simulation.sjon in the pngine repository.

How it works

Section titled “How it works”

The document is SJON, the S-expression format pngine compiles: each form is one WebGPU resource or operation, and the shader text inside (shader-module …) is plain WGSL, WebGPU’s shading language. The ping-pong structure is that of Game of Life, doubled: two :pool 2 buffers instead of one.

Two fields, four bindings

Section titled “Two fields, four bindings”

velocityBuffer is 96 × 96 × 2 floats and densityBuffer 96 × 96 floats, both [storage] with :pool 2. The step bind group therefore has five entries: the uniforms, then velIn / velOut and denIn / denOut with :ping-pong 0 / :ping-pong 1 on each pair. One :bind-groups-pool-offsets [0] on the pass, advanced by the runtime every frame, flips both pairs together. (init :name initFluid …) zeroes the velocity once with a compute shader; the density relies on WebGPU’s guarantee that a new buffer starts zeroed.

Semi-Lagrangian advection

Section titled “Semi-Lagrangian advection”

For each cell, stepShader reads its velocity, steps backwards by vel · DT (DT = 0.5), and asks what the fields were at that point: sampleVelBilinear and sampleDenBilinear clamp the coordinates to the grid, take the four surrounding cells and mix them by the fractional part. Pulling values from upstream instead of pushing them downstream is what makes the scheme unconditionally stable: nothing can overshoot.

Two sources then add force and dye inside discs whose centres orbit at different rates: one pushes along a rotating direction with strength 8, the other with strength 6, both falling off linearly to the disc’s edge, and each deposits density in the same disc. Cells within two of the boundary have their velocity halved, and both fields decay by 0.995 per step; the density is clamped to [0, 3].

(compute-pass :name stepPass … (dispatch :workgroups [6 6 1])) covers the grid with 6 × 16 = 96 threads per axis.

The display

Section titled “The display”

renderShader binds the uniforms plus the density and velocity variants just written (:ping-pong 1). Per pixel it finds the cell, takes the velocity’s angle with atan2 as a hue and its speed as value through hsv2rgb, weights that by 0.6, and adds the density as a grey brightness. Direction shows as colour, speed as intensity, dye as white.

In the specifications

Section titled “In the specifications”
What the sample uses WebGPU WGSL
Two ping-pong storage buffers in one bind group bind group creation, GPUBindGroupEntry, GPUBufferUsage.STORAGE storage address space, access modes, runtime-sized arrays
Compute pass over the grid compute passes, dispatchWorkgroups() @workgroup_size, global_invocation_id
Bilinear sampling by hand floor, fract, clamp, min, mix, i32() / u32() conversions
Fullscreen display draw() vertex_index, atan2, length
Sources driven by time writeBuffer() uniform address space, sin / cos
Section titled “Related”