Wave simulation

Download PNG

Ripples on a pond. Three oscillators wander over a 128×128 height field and push it up and down; each frame a compute shader integrates the wave equation one step, and a fullscreen triangle lights the surface from the slope of the height, so crests catch a specular highlight and troughs go dark blue. Unlike diffusion, the wave equation is second-order in time, so every cell keeps two values: its current height and its previous one. The poster is the first, still frame; press Play.

; Wave simulation. The 2D wave equation on a 128×128 height field: three
; moving oscillators excite it, a compute pass integrates the finite-difference
; Laplacian with damping each frame, and a fullscreen triangle shades the
; height with gradient-derived normals (diffuse + specular) in water colours.
; The field is one `:pool 2` ping-pong buffer of `array<vec2f>` (.x current
; height, .y previous height); a compute `(init …)` zeroes it once, the step
; bind group has three entries (uniforms + wave in/out with `:ping-pong`),
; passes select the variant with `:bind-groups-pool-offsets`, the step runs
; `(dispatch :workgroups [8 8 1])` over a 16×16 workgroup, and `fs` reads the storage buffer
; directly.
;
; GRID_SIZE sizes the buffer; the shaders do not read it. `(define …)` values
; are expression-only and never reach WGSL text, so each shader restates the
; grid extent as a local `const SIZE: u32 = 128u;` / `const SIZE: f32 = 128.0;`.
; Both spellings must be kept in step by hand: change one, change the other.

(define :name GRID_SIZE :value 128)

; The ping-pong wave field: `:pool 2` allocates two variants.
(buffer :name waveBuffer :size (* GRID_SIZE GRID_SIZE 2 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 wave field once.
(init :name initWave :buffer waveBuffer :module initShader :workgroups [(ceil (/ (* GRID_SIZE GRID_SIZE) 64))])

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

@binding(0) @group(0) var<storage, read_write> wave: Wave;

const SIZE: u32 = 128u;

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

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

struct Wave { data: array<vec2f> }

@group(0) @binding(1) var<storage, read> waveIn: Wave;
@group(0) @binding(2) var<storage, read_write> waveOut: Wave;

const SIZE: u32 = 128u;
const SIZEF: f32 = 128.0;
const C: f32 = 0.3;
const DAMPING: f32 = 0.998;

fn getHeight(x: i32, y: i32) -> f32 {
  if (x < 0 || x >= i32(SIZE) || y < 0 || y >= i32(SIZE)) {
    return 0.0;
  }
  return waveIn.data[u32(y) * SIZE + u32(x)].x;
}

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

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

  let state = waveIn.data[idx];
  let h = state.x;       // current
  let hPrev = state.y;   // previous

  // Wave equation with Laplacian
  let laplacian = getHeight(x-1, y) + getHeight(x+1, y) +
                  getHeight(x, y-1) + getHeight(x, y+1) - 4.0 * h;

  var hNew = 2.0 * h - hPrev + C * C * laplacian;

  // Add wave sources
  let fx = f32(id.x);
  let fy = f32(id.y);

  // Source 1: oscillating position
  let cx1 = SIZEF * 0.5 + sin(u.time * 1.2) * SIZEF * 0.3;
  let cy1 = SIZEF * 0.5 + cos(u.time * 0.9) * SIZEF * 0.3;
  let d1 = length(vec2f(fx - cx1, fy - cy1));
  if (d1 < 2.5) {
    hNew += sin(u.time * 15.0) * 0.3;
  }

  // Source 2: different frequency
  let cx2 = SIZEF * 0.3 + sin(u.time * 0.8 + 2.0) * SIZEF * 0.2;
  let cy2 = SIZEF * 0.7 + cos(u.time * 1.1 + 1.0) * SIZEF * 0.2;
  let d2 = length(vec2f(fx - cx2, fy - cy2));
  if (d2 < 2.5) {
    hNew += sin(u.time * 12.0 + 1.5) * 0.25;
  }

  // Source 3: slower, larger waves
  let cx3 = SIZEF * 0.7 + sin(u.time * 0.5) * SIZEF * 0.15;
  let cy3 = SIZEF * 0.3 + cos(u.time * 0.7) * SIZEF * 0.15;
  let d3 = length(vec2f(fx - cx3, fy - cy3));
  if (d3 < 3.0) {
    hNew += sin(u.time * 8.0) * 0.2;
  }

  // Damping and clamping
  hNew *= DAMPING;
  hNew = clamp(hNew, -1.0, 1.0);

  // Store: new becomes current, current becomes previous
  waveOut.data[idx] = vec2f(hNew, h);
}
""")

(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 waveBuffer :ping-pong 0)
  (entry :binding 2 :buffer waveBuffer :ping-pong 1))

(compute-pass :name stepPass
  :pipeline stepPipeline
  :bind-groups [stepBindGroup]
  :bind-groups-pool-offsets [0]
  (dispatch :workgroups [8 8 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 Wave { data: array<vec2f> }
@group(0) @binding(1) var<storage, read> wave: Wave;

const SIZE: f32 = 128.0;
const SIZEI: u32 = 128u;

@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 getH(x: u32, y: u32) -> f32 {
  if (x >= SIZEI || y >= SIZEI) { return 0.0; }
  return wave.data[y * SIZEI + x].x;
}

@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 >= SIZEI || cellY >= SIZEI) {
    return vec4f(0.0, 0.05, 0.1, 1.0);
  }

  let h = getH(cellX, cellY);

  // Compute gradient for lighting
  let hL = getH(max(0u, cellX - 1u), cellY);
  let hR = getH(min(SIZEI - 1u, cellX + 1u), cellY);
  let hD = getH(cellX, max(0u, cellY - 1u));
  let hU = getH(cellX, min(SIZEI - 1u, cellY + 1u));

  let gradX = (hR - hL) * 0.5;
  let gradY = (hU - hD) * 0.5;

  // Simple normal and lighting
  let normal = normalize(vec3f(-gradX * 3.0, -gradY * 3.0, 1.0));
  let lightDir = normalize(vec3f(0.3, 0.4, 1.0));
  let diff = max(dot(normal, lightDir), 0.2);

  // Specular
  let viewDir = vec3f(0.0, 0.0, 1.0);
  let halfDir = normalize(lightDir + viewDir);
  let spec = pow(max(dot(normal, halfDir), 0.0), 32.0);

  // Water color based on height
  let waterDeep = vec3f(0.0, 0.1, 0.25);
  let waterMid = vec3f(0.1, 0.35, 0.5);
  let waterLight = vec3f(0.3, 0.6, 0.8);
  let foam = vec3f(0.8, 0.9, 1.0);

  let t = h * 0.5 + 0.5;
  var color = mix(waterDeep, waterMid, smoothstep(0.3, 0.5, t));
  color = mix(color, waterLight, smoothstep(0.5, 0.7, t));
  color = mix(color, foam, smoothstep(0.75, 0.95, t));

  color = color * diff + vec3f(1.0, 0.95, 0.9) * spec * 0.4;

  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 waveBuffer :ping-pong 1))

(render-pass :name drawPass
  (color-attachment :view context-current-texture :clear-value [0 0.05 0.1 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 [initWave]
  :perform [writeUniforms stepPass drawPass])

examples/samples/22_wave_simulation.sjon in the pngine repository.

.sample-source[data-folded] .expressive-code { max-block-size: none !important; mask-image: none !important; } .sample-source[data-folded] .sample-source-more { display: none; }

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 structure is that of Game of Life: one :pool 2 buffer, a step pass with :ping-pong bindings, a fullscreen display reading the variant just written.

Two heights per cell

(buffer :name waveBuffer :size (* GRID_SIZE GRID_SIZE 2 4) :usage [storage] :pool 2) holds an array<vec2f>: .x the current height, .y the previous. That is what lets a single buffer carry a second-order scheme, and it is why the step ends with waveOut.data[idx] = vec2f(hNew, h): the new height becomes current, the current becomes previous. (init :name initWave …) zeroes both once. As on the other simulation pages, (define :name GRID_SIZE …) sizes the buffer only, and each shader restates the extent as const SIZE, once as u32 and once as f32, which have to be kept in step by hand.

The step

Per cell, with h and hPrev from waveIn and neighbours through getHeight (0 outside the grid: a fixed, reflecting boundary):

laplacian = left + right + up + down - 4h;
hNew = 2h - hPrev + C² · laplacian, the leapfrog form of the wave equation, with C = 0.3 well under the stability limit;
three sources add sin(u.time · f) · amplitude inside small discs whose centres orbit at different rates and radii, at frequencies 15, 12 and 8;
hNew *= DAMPING (0.998) and clamp to [-1, 1].

(compute-pass :name stepPass … :bind-groups-pool-offsets [0] (dispatch :workgroups [8 8 1])) runs 8 × 8 workgroups of 16 × 16 threads; the runtime advances the pool offset every frame so the roles of the two variants swap.

Shading from the slope

renderShader reads the height at the pixel’s cell and its four neighbours, forms a central-difference gradient, and builds a normal normalize(vec3f(-3 gradX, -3 gradY, 1)), the factor 3 exaggerating the slope. From there it is Blinn-Phong: max(dot(N, L), 0.2) diffuse and a pow(…, 32) half-vector specular under a fixed light. The base colour is a three-stop water ramp by height (deep, mid, light) with a foam colour mixed in near the top of the range, so crests turn pale.

In the specifications

What the sample uses	WebGPU	WGSL
Storage buffer of `vec2f`, ping-ponged	`GPUBufferUsage.STORAGE`, bind group creation	`storage` address space, runtime-sized arrays, vector types
Compute pass over the grid	compute passes, `dispatchWorkgroups()`	`@workgroup_size`, `global_invocation_id`
Fullscreen display	`draw()`	`vertex_index`, `position`
Time for the oscillators	`writeBuffer()`	`uniform` address space, `sin` / `cos`
Integration and shading		`clamp`, `normalize`, `dot`, `pow`, `smoothstep`, `mix`, `max` / `min`

Diffusion is the first-order cousin (one value per cell); Game of Life explains the ping-pong forms in detail; Simple lighting is the same Blinn-Phong on geometry.
Forms: (buffer …) (:pool), (bind-group …) (:ping-pong), (compute-pass …), (define …).