Lattice Boltzmann Fluid in WebGL
Run a real fluid simulation entirely on the GPU. Using a pair of floating-point textures and a fragment shader, each pixel computes one lattice node's distribution functions — then streams, collides, and visualises the fluid velocity field in real time.
1. LBM Theory — D2Q9
The Lattice Boltzmann Method doesn't solve the Navier-Stokes equations directly. Instead it simulates a gas of particles on a grid where each node holds 9 distribution functions f_i (the D2Q9 model), representing particle populations moving in 9 discrete velocity directions.
The 9 velocity vectors (plus rest) and their weights:
Direction: e_i weight w_i rest: (0,0) 4/9 E,N,W,S: (±1,0),(0,±1) 1/9 NE,NW,SW,SE: (±1,±1) 1/36
Each step has two phases: streaming (populations advect along their velocity direction) and collision (populations relax toward a local equilibrium, controlled by viscosity):
// BGK collision (single relaxation time) f_eq_i = w_i * ρ * (1 + 3(e_i·u) + 4.5(e_i·u)² − 1.5(u·u)) f_i* = f_i + (f_eq_i − f_i) / τ // collision result // macroscopic quantities: ρ = Σ f_i // density ρu = Σ f_i * e_i // momentum ν = (τ − 0.5) / 3 // kinematic viscosity
2. Texture Layout
D2Q9 needs 9 floats per node. We use two RGBA_FLOAT textures: one holds f0..f3 packed as RGBA (rest + E + N + W), and one holds f4..f7 (S + NE + NW + SW) plus f8 in the alpha channel, with the obstacle mask packed in separately.
// Texture A: RGBA float (width × height) // R = f0 (rest) // G = f1 (E) // B = f2 (N) // A = f3 (W) // Texture B: RGBA float // R = f4 (S) // G = f5 (NE) // B = f6 (NW) // A = f7 (SW) // Texture C: f8 (SE) + obstacle mask // R = f8 (SE) // G = obstacle (0 or 1) // B, A unused
OES_texture_float and WEBGL_color_buffer_float extensions. WebGL2 supports float textures natively. Use WebGL2 for new projects.3. WebGL Setup
- Create context and extensions
const canvas = document.querySelector('canvas'); const gl = canvas.getContext('webgl2'); if (!gl) throw new Error('WebGL2 required'); // Float texture support (WebGL2 has it by default) // But MRT (multiple render targets) needs: // gl.getExtension('EXT_color_buffer_float'); - Create ping-pong framebuffers
function createFBOTexture(gl, w, h) { const tex = gl.createTexture(); gl.bindTexture(gl.TEXTURE_2D, tex); gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA32F, w, h, 0, gl.RGBA, gl.FLOAT, null); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); const fbo = gl.createFramebuffer(); gl.bindFramebuffer(gl.FRAMEBUFFER, fbo); gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0); return { tex, fbo }; } // Create two sets of A/B/C textures (current and next) const stateA = [createFBOTexture(gl, W, H), createFBOTexture(gl, W, H)]; let readIdx = 0, writeIdx = 1; - Full-screen quad for compute passes
// Vertex shader: just a full-screen triangle const vsSource = `#version 300 es in vec2 aPos; out vec2 vUV; void main() { vUV = aPos * 0.5 + 0.5; gl_Position = vec4(aPos, 0.0, 1.0); }`; // Two triangles covering NDC: const positions = new Float32Array([ -1,-1, 1,-1, -1,1, -1, 1, 1,-1, 1,1 ]);
4. Streaming Shader
The streaming pass samples from neighbours. For direction i with velocity (cx, cy), f_i at position (x,y) comes from position (x−cx, y−cy) in the previous step:
#version 300 es
precision highp float;
uniform sampler2D uTexA; // f0..f3
uniform sampler2D uTexB; // f4..f7
uniform sampler2D uTexC; // f8 + obstacle
uniform vec2 uInvRes;
in vec2 vUV;
out vec4 fragColor; // updated texA
vec2 e[9] = vec2[](
vec2(0,0), vec2(1,0), vec2(0,1), vec2(-1,0),
vec2(0,-1), vec2(1,1), vec2(-1,1), vec2(-1,-1), vec2(1,-1)
);
float sampleF(int i, vec2 uv) {
vec2 offset = -e[i] * uInvRes; // reverse: pull from upstream
vec4 a = texture(uTexA, uv + offset);
vec4 b = texture(uTexB, uv + offset);
float c = texture(uTexC, uv + offset).r;
// extract channel for direction i
if (i == 0) return a.r;
if (i == 1) return a.g;
if (i == 2) return a.b;
if (i == 3) return a.a;
if (i == 4) return b.r;
if (i == 5) return b.g;
if (i == 6) return b.b;
if (i == 7) return b.a;
return c; // i == 8
}
void main() {
// Stream: f0..f3 into texA
float f0 = sampleF(0, vUV);
float f1 = sampleF(1, vUV);
float f2 = sampleF(2, vUV);
float f3 = sampleF(3, vUV);
fragColor = vec4(f0, f1, f2, f3);
}5. Collision Shader
After streaming, the collision pass computes macroscopic ρ and u, then relaxes toward equilibrium:
#version 300 es
precision highp float;
uniform sampler2D uTexA; // streamed f0..f3
uniform sampler2D uTexB; // streamed f4..f7
uniform sampler2D uTexC; // streamed f8
uniform float uTau; // relaxation time (1.5–2.0 for stability)
in vec2 vUV;
out vec4 fragColor;
const float W[9] = float[](4./9., 1./9., 1./9., 1./9., 1./9.,
1./36., 1./36., 1./36., 1./36.);
const vec2 E[9] = vec2[](
vec2(0,0),vec2(1,0),vec2(0,1),vec2(-1,0),vec2(0,-1),
vec2(1,1),vec2(-1,1),vec2(-1,-1),vec2(1,-1));
float feq(int i, float rho, vec2 u) {
float eu = dot(E[i], u);
float uu = dot(u, u);
return W[i] * rho * (1.0 + 3.0*eu + 4.5*eu*eu - 1.5*uu);
}
void main() {
vec4 a = texture(uTexA, vUV);
vec4 b = texture(uTexB, vUV);
float f8 = texture(uTexC, vUV).r;
float f[9];
f[0]=a.r; f[1]=a.g; f[2]=a.b; f[3]=a.a;
f[4]=b.r; f[5]=b.g; f[6]=b.b; f[7]=b.a;
f[8]=f8;
// Macroscopic quantities
float rho = 0.0;
vec2 u = vec2(0.0);
for (int i = 0; i < 9; i++) {
rho += f[i];
u += f[i] * E[i];
}
u /= rho;
// BGK collision
float inv_tau = 1.0 / uTau;
for (int i = 0; i < 9; i++) {
f[i] += (feq(i, rho, u) - f[i]) * inv_tau;
}
fragColor = vec4(f[0], f[1], f[2], f[3]);
// write texB and texC similarly via MRT layout
}6. Obstacle Bounce-Back
Obstacles are stored as a binary mask (1 = solid) in the green channel of texC. When a node is solid, instead of normal streaming+collision, apply no-slip bounce-back: reflect each distribution function to its opposite direction:
// In streaming shader, after sampling streamed f:
float isObstacle = texture(uTexC, vUV).g;
if (isObstacle > 0.5) {
// Bounce-back: swap opposite directions
// f1(E) ↔ f3(W), f2(N) ↔ f4(S), f5(NE) ↔ f7(SW), f6(NW) ↔ f8(SE)
float tmp;
tmp = f[1]; f[1] = f[3]; f[3] = tmp;
tmp = f[2]; f[2] = f[4]; f[4] = tmp;
tmp = f[5]; f[5] = f[7]; f[7] = tmp;
tmp = f[6]; f[6] = f[8]; f[8] = tmp;
}Draw obstacles interactively: in a JavaScript mousemove handler, render filled circles to the obstacle texture. The bounce-back correctly enforces no-slip boundary conditions at any obstacle shape.
7. Velocity Visualisation
The final display pass reads the current fluid state, computes ρ and u, then maps velocity to colour:
// Visualisation fragment shader
void main() {
vec4 a = texture(uTexA, vUV);
vec4 b = texture(uTexB, vUV);
float f8 = texture(uTexC, vUV).r;
float obs = texture(uTexC, vUV).g;
float f[9];
f[0]=a.r;f[1]=a.g;f[2]=a.b;f[3]=a.a;
f[4]=b.r;f[5]=b.g;f[6]=b.b;f[7]=b.a;f[8]=f8;
vec2 u = vec2(0.0);
float rho = 0.0;
for (int i=0; i<9; i++) { rho+=f[i]; u+=f[i]*E[i]; }
u /= rho;
float speed = length(u);
// HSV: hue = velocity direction, saturation+value = speed
float hue = atan(u.y, u.x) / (2.0 * PI) + 0.5;
vec3 col = hsv2rgb(vec3(hue, 0.85, clamp(speed * 20.0, 0.0, 1.0)));
if (obs > 0.5) col = vec3(0.2); // dark grey obstacle
outColor = vec4(col, 1.0);
}