Lattice-Boltzmann Fluid in 200 Lines

D2Q9 = 2 dimensions, 9 velocity vectors. Each lattice cell stores 9 distribution values f[0..8] — the probability density of particles moving in each direction:

// D2Q9 velocity directions: [ex, ey] for i = 0..8
//   6  2  5
//   3  0  1    (0 = rest, 1-4 = cardinal, 5-8 = diagonal)
//   7  4  8
const EX = [ 0,  1,  0, -1,  0,  1, -1, -1,  1];
const EY = [ 0,  0,  1,  0, -1,  1,  1, -1, -1];
const W  = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; // weights
const OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]; // opposite direction index

const NX = 200, NY = 100;
const N  = NX * NY;

// 9 distributions per cell, stored contiguously: f[i*N + (y*NX + x)]
const f    = new Float64Array(9 * N);
const fNew = new Float64Array(9 * N);
const wall = new Uint8Array(N);   // 1 = solid obstacle

function at(x, y) { return y * NX + x; }

Initialise every cell to the equilibrium distribution for a uniform horizontal flow at speed u0:

const u0 = 0.1; // initial flow speed (in lattice units, keep < 0.3)

function feq(i, rho, ux, uy) {
  const eu = EX[i] * ux + EY[i] * uy;
  const u2 = ux*ux + uy*uy;
  return W[i] * rho * (1 + 3*eu + 4.5*eu*eu - 1.5*u2);
}

function initFlow() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      for (let i = 0; i < 9; i++) {
        f[i*N + cell] = feq(i, 1.0, u0, 0.0);
      }
    }
  }
}
initFlow();

Stream: each distribution f[i] moves to the neighbour cell in direction i. This is just an array copy pattern:

function stream() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      for (let i = 0; i < 9; i++) {
        // Source: where did this distribution come FROM?
        const sx = x - EX[i], sy = y - EY[i];
        if (sx < 0 || sx >= NX || sy < 0 || sy >= NY) {
          fNew[i*N + cell] = f[i*N + cell]; // edge: keep
        } else {
          fNew[i*N + cell] = f[i*N + at(sx, sy)];
        }
      }
    }
  }
  f.set(fNew);
}

The BGK collision operator relaxes f towards its equilibrium feq. The relaxation time tau controls viscosity:

// Kinematic viscosity: nu = (tau - 0.5) / 3
// tau = 1.0 → nu = 1/6 (moderate viscosity)
// tau = 0.6 → nu = 1/30 (lower viscosity, Re ≈ u0 * L / nu)
const tau = 0.6;
const omega = 1.0 / tau; // = 1/tau

function collide() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      if (wall[cell]) continue; // skip solid cells

      // Compute macroscopic density and velocity
      let rho = 0, ux = 0, uy = 0;
      for (let i = 0; i < 9; i++) {
        const fi = f[i*N + cell];
        rho += fi;
        ux  += EX[i] * fi;
        uy  += EY[i] * fi;
      }
      ux /= rho; uy /= rho;

      // Relax towards equilibrium
      for (let i = 0; i < 9; i++) {
        const eq = feq(i, rho, ux, uy);
        f[i*N + cell] += omega * (eq - f[i*N + cell]);
      }
    }
  }
}

For solid walls: particles that stream into a wall cell are reflected back the way they came (no-slip condition):

// Add a circular obstacle
function addCircle(cx, cy, radius) {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if ((x-cx)**2 + (y-cy)**2 < radius**2) wall[at(x,y)] = 1;
    }
  }
}
addCircle(NX * 0.3, NY * 0.5, NY * 0.1);

// Apply bounce-back AFTER streaming
function bounceBack() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      if (!wall[cell]) continue;
      for (let i = 0; i < 9; i++) {
        // Swap with opposite direction
        const opp = OPP[i];
        const tmp = f[i*N + cell];
        f[i*N + cell] = f[opp*N + cell];
        f[opp*N + cell] = tmp;
      }
    }
  }
}

Reset the left column to inlet flow and zero-gradient the right column each step:

function applyBoundaries() {
  // Left inlet: fixed velocity u0
  for (let y = 0; y < NY; y++) {
    const cell = at(0, y);
    for (let i = 0; i < 9; i++) f[i*N + cell] = feq(i, 1.0, u0, 0.0);
  }
  // Right outlet: copy from x = NX-2 (zero normal-gradient)
  for (let y = 0; y < NY; y++) {
    const c = at(NX-1, y), c2 = at(NX-2, y);
    for (let i = 0; i < 9; i++) f[i*N + c] = f[i*N + c2];
  }
}

// Main simulation loop
function step() {
  stream();
  bounceBack();
  collide();
  applyBoundaries();
}

Curl (∂uy/∂x − ∂ux/∂y) highlights vortices beautifully. Use ImageData to write pixels directly — much faster than fillRect:

const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
canvas.width = NX; canvas.height = NY;
const img = ctx.createImageData(NX, NY);

function getVelocity(x, y) {
  const cell = at(x, y);
  let ux = 0, uy = 0, rho = 0;
  for (let i = 0; i < 9; i++) {
    const fi = f[i*N + cell];
    rho += fi; ux += EX[i]*fi; uy += EY[i]*fi;
  }
  return { ux: ux/rho, uy: uy/rho };
}

function render() {
  for (let y = 1; y < NY-1; y++) {
    for (let x = 1; x < NX-1; x++) {
      const p = (y * NX + x) * 4;

      if (wall[at(x, y)]) {
        img.data[p] = img.data[p+1] = img.data[p+2] = 60;
        img.data[p+3] = 255;
        continue;
      }

      // Finite difference curl
      const r = getVelocity(x+1, y), l = getVelocity(x-1, y);
      const t = getVelocity(x, y+1), b = getVelocity(x, y-1);
      const curl = (r.uy - l.uy) - (t.ux - b.ux);

      // Map curl to colour: negative = blue, positive = red
      const v = Math.max(-0.1, Math.min(0.1, curl)) / 0.1; // [-1, 1]
      img.data[p]   = v > 0 ? Math.floor(255 * v) : 0;
      img.data[p+1] = 0;
      img.data[p+2] = v < 0 ? Math.floor(-255 * v) : 0;
      img.data[p+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}

// Animate: run multiple physics steps per frame for speed
let running = true;
(function loop() {
  if (!running) return;
  for (let s = 0; s < 10; s++) step(); // 10 sub-steps per frame
  render();
  requestAnimationFrame(loop);
})();