@vshaders/fluid

Stable Fluids simulation kernels as pure WGSL modules for vgpu shaders: semi-Lagrangian advection, pressure projection (divergence, Jacobi relaxation, gradient subtraction), vorticity confinement, and Gaussian splats. Each kernel is one pass of the classic Stam splitting; entry shaders own the bindings and pass textures in:

import { advectBacktrace, bilerp2, bilerp4 } from "@vshaders/fluid/advect";

struct Uniforms {
  resolution: vec2f,
  dt: f32,
  dissipation: f32,
}

@group(0) @binding(0) var<uniform> uniforms: Uniforms;
@group(0) @binding(1) var srcTex: texture_2d<f32>;
@group(0) @binding(2) var velocityTex: texture_2d<f32>;
@group(0) @binding(3) var linearSampler: sampler;

@fragment fn fs_main(@location(0) uv: vec2f) -> @location(0) vec4f {
  let texel = 1.0 / max(uniforms.resolution, vec2f(1.0));
  let velocity = bilerp2(velocityTex, linearSampler, uv);
  let back = advectBacktrace(uv, velocity, uniforms.dt, texel);
  return bilerp4(srcTex, linearSampler, back) * uniforms.dissipation;
}

Every module is pure WGSL — functions only, no bindings, no entry points. Functions take texture_2d<f32> and sampler parameters where they sample fields; the entry shader declares the bindings and passes them through, so vgpu's resolver composes and prunes as usual.

Conventions

The kernels agree on one grid model; mix them freely as long as your fields do too.

• Fields. velocityTex stores 2D velocity in .xy (an rg16float target works). pressureTex, divergenceTex, and curlTex store scalars in .x (r16float). Dye is any four-component field.
• Units. uv spans [0, 1] across a field; texel is one grid cell in uv, i.e. 1.0 / resolution of the field being processed. Velocity is stored in grid texels per unit of dt. Grid spacing is one texel, so all central differences carry the 0.5 factor and divergence, pressure, and gradients share per-texel units.
• Boundaries. Integer-coordinate kernels (*At) read neighborhoods with textureLoad at coordinates clamped to the texture extent — clamp-to-edge, approximating free-slip walls, with no sampler and no filterability requirement. Sampling kernels expect a linear-filtering, clamp-to-edge sampler for the same boundary behavior. The projection kernels also come in *BoundedAt variants with reflective container walls; see @vshaders/fluid/project.
• Coordinates. p: vec2i is the integer texel coordinate of the cell being written — vec2i(position.xy) in a fullscreen fragment pass. Out-of-range p clamps like every other read.

@vshaders/fluid/advect

Semi-Lagrangian advection: to move a field through the flow, look backward along the velocity and take what was there. Linear interpolation at the backtraced point is what makes the scheme unconditionally stable.

• bilerp2(tex: texture_2d<f32>, samp: sampler, uv: vec2f) -> vec2f — bilinear sample of a two-component field (velocity) at uv. samp must filter linearly for true bilinear interpolation.
• bilerp4(tex: texture_2d<f32>, samp: sampler, uv: vec2f) -> vec4f — the four-component variant, for dye.
• advectBacktrace(uv: vec2f, velocity: vec2f, dt: f32, texel: vec2f) -> vec2f — where the quantity now at uv was one step ago: uv - velocity * dt * texel. velocity is the flow at uv (from bilerp2 on the velocity field); texel converts it from texels per unit dt into uv displacement. Non-positive texel components contribute no displacement (a degenerate axis advects in place) instead of reversing the trace. Sample the advected field at the result and scale by your dissipation factor.

@vshaders/fluid/project

The projection that makes the flow incompressible: measure divergence, relax the pressure Poisson equation, subtract the pressure gradient. All three run on the velocity grid and are numerically consistent with each other.

• divergenceAt(velocityTex: texture_2d<f32>, p: vec2i) -> f32 — central-difference divergence of the velocity field at p. Write it once per projection into a scalar target.
• jacobiPressureAt(pressureTex: texture_2d<f32>, divergenceTex: texture_2d<f32>, p: vec2i) -> f32 — one Jacobi relaxation step of laplacian(pressure) = divergence: the next pressure at p from the previous iteration (pressureTex) and the fixed divergence. Iterate by ping-ponging the pressure texture (20 iterations is a good default for pointer-splat-scale flows; low-frequency divergence converges much more slowly).
• pressureGradientAt(pressureTex: texture_2d<f32>, p: vec2i) -> vec2f — central-difference gradient of the relaxed pressure.
• subtractPressureGradientAt(velocityTex: texture_2d<f32>, pressureTex: texture_2d<f32>, p: vec2i) -> vec2f — the projected velocity: velocity - pressureGradient at p, approximately divergence-free. Write it back to the velocity field.

Boundary conditions

The plain kernels read out-of-range neighbors clamped to the edge (clamp-to-edge), under which nothing stops flow at the domain edge: fluid drifts out and the walls are effectively open. Each kernel has a *BoundedAt variant with the same signature that treats the domain edge as a reflective (no-through-flow) container wall instead, using the ghost-cell boundary conditions of Harris, GPU Gems ch. 38:

• Velocity (divergenceBoundedAt): at a boundary texel the off-grid neighbor's velocity is the negation of the center's, so the wall-face velocity — the average of the two — is zero. Flow into a wall registers as compression in the divergence, and the pressure solve pushes back: incoming fluid bounces.
• Pressure (jacobiPressureBoundedAt, pressureGradientBoundedAt, subtractPressureGradientBoundedAt): the off-grid pressure neighbor equals the center — pure Neumann, dp/dn = 0 at walls — so the projection never accelerates flow through a wall. On a same-size grid this coincides numerically with the clamped free reads (the clamped edge read is the center); the bounded variants make the wall condition explicit and keep a contained projection reading as one family.

Use one family per projection: divergenceBoundedAt + jacobiPressureBoundedAt + subtractPressureGradientBoundedAt for a contained domain, the plain kernels for a free one. Since the signatures match, a contained uniform can select between them per pass (a select on a uniform keeps control flow uniform).

• divergenceBoundedAt(velocityTex: texture_2d<f32>, p: vec2i) -> f32 — divergence with reflective walls. The wall sits half a texel outside the outermost texel ring; a one-texel grid axis degenerates gracefully (both ghosts cancel, zero difference along that axis).
• jacobiPressureBoundedAt(pressureTex: texture_2d<f32>, divergenceTex: texture_2d<f32>, p: vec2i) -> f32 — one Jacobi step with pure-Neumann walls.
• pressureGradientBoundedAt(pressureTex: texture_2d<f32>, p: vec2i) -> vec2f — pressure gradient with pure-Neumann walls: zero normal derivative at the boundary.
• subtractPressureGradientBoundedAt(velocityTex: texture_2d<f32>, pressureTex: texture_2d<f32>, p: vec2i) -> vec2f — the contained projection step.

Advection needs no bounded variant. Semi-Lagrangian advection with a clamp-to-edge sampler backtraces past a wall to the edge value, and the bounded projection keeps the wall-normal flow near zero, so backtraces rarely leave the domain in the first place. Verified empirically: with the bounded projection kernels and unchanged advection, a stroke driven straight at a wall mushrooms back along it instead of exiting, and the domain retains its dye (wall-column outflow collapses to ~0 after impact, where the free kernels leave it growing).

@vshaders/fluid/vorticity

Grid advection dissipates small vortices; vorticity confinement (Fedkiw et al. 2001) measures the curl the grid still has and steers flow back around it.

• curlAt(velocityTex: texture_2d<f32>, p: vec2i) -> f32 — scalar (out-of-plane) curl of the velocity field at p by central differences. Write it into a scalar target before the confinement pass.
• vorticityForceAt(curlTex: texture_2d<f32>, p: vec2i, strength: f32) -> vec2f — the confinement force at p from the curl field: strength * curl along the rotated unit gradient of |curl|. Add force * dt to velocity. Returns zero where the |curl| gradient vanishes (no vortex structure to sharpen) instead of normalizing a zero vector, and for non-positive strength. Typical strength is 5–30 on a 128-grid; more is stormier.

@vshaders/fluid/splat

• gaussianSplat(p: vec2f, center: vec2f, radius: f32) -> f32 — the Gaussian falloff exp(-|p - center|^2 / radius^2) used to inject force and dye: 1 at the center, 1/e at distance radius, effectively zero past ~3 radii. p, center, and radius must share one coordinate space (correct uv for aspect ratio before calling, or work in pixels). A non-positive radius deposits nothing instead of dividing by zero. Scale the falloff by your force or dye color and add it to the field being splatted.

A projection frame

The kernels compose into the classic per-tick pass graph — advect, (splat, confine,) project — as fullscreen fragment passes over ping-pong float targets, one frame() per tick:

import { effect, frame, init, pingPong, sampler, target } from "vgpu";

const gpu = await init();
const velocity = pingPong(gpu, 128, 128, { format: "rg16float" });
const pressure = pingPong(gpu, 128, 128, { format: "r16float" });
const divergence = target(gpu, { size: [128, 128], format: "r16float" });
const linear = sampler(gpu, {
  minFilter: "linear", magFilter: "linear",
  addressModeU: "clamp-to-edge", addressModeV: "clamp-to-edge",
});

// advect, div, jacobi, subtract = effect(gpu, entryShader) for entries
// importing the kernels above.
frame(gpu, (f) => {
  advect.set({ srcTex: velocity.read, velocityTex: velocity.read, linearSampler: linear });
  f.pass(velocity.write, advect);
  velocity.swap();

  div.set({ velocityTex: velocity.read });
  f.pass(divergence, div);

  for (let i = 0; i < 20; i++) {
    jacobi.set({ pressureTex: pressure.read, divergenceTex: divergence });
    f.pass(pressure.write, jacobi);
    pressure.swap();
  }

  subtract.set({ velocityTex: velocity.read, pressureTex: pressure.read });
  f.pass(velocity.write, subtract);
  velocity.swap();
});

Texture rebinds between passes of one frame take effect per pass; JS-value uniform writes do not (they are frame-global), so per-pass uniform values need one effect instance each.

Provenance

The method is Jos Stam's "Stable Fluids" (SIGGRAPH 1999) and its game-oriented restatement "Real-Time Fluid Dynamics for Games" (GDC 2003): semi-Lagrangian advection plus pressure projection via Jacobi relaxation. Vorticity confinement is from Fedkiw, Stam, and Jensen, "Visual Simulation of Smoke" (SIGGRAPH 2001). The texture-based GPU formulation (fields as textures, kernels as fragment passes, central differences at one-texel spacing) follows Mark Harris, "Fast Fluid Dynamics Simulation on the GPU", GPU Gems ch. 38 (2004). The Gaussian splat is the standard injection used throughout that literature. All WGSL in this package is an original implementation from the published equations — no shader code was transcribed from any existing fluid demo — with edge-case guards (clamped neighborhoods, zero-gradient confinement, non-positive radius, strength, and texel sizes) in the argument ranges the math leaves undefined.

Verifying

npx vgpu check path/to/your-entry-shader.wgsl

resolves the import graph, validates the composed shader, and prints its reflection.

License

MIT