SPH Fluid Simulation from Scratch in 200 Lines

// Smoothing radius const H = 30; // px — all neighbours must be within this distance const H2 = H * H; // Poly6 kernel — for density estimation (smooth, even at r=0) function Wpoly6(r2) { if (r2 >= H2) return 0; const x = H2 - r2; return 4 / (Math.PI * H2 * H2 * H2 * H2) * x * x * x; // exact prefactor: 4/(pi * h^8), but we bake into REST_DENSITY choice } // Spiky kernel gradient — for PRESSURE forces (non-zero at r=0, resists clumping) function gradWspiky(rx, ry, r) { if (r <= 0 || r >= H) return { x: 0, y: 0 }; const coeff = -30 / (Math.PI * H2 * H2 * H2) * (H - r) * (H - r) / r; return { x: coeff * rx, y: coeff * ry }; } // Viscosity kernel Laplacian — for VISCOSITY forces (positive, diffuses velocity) function lapWviscosity(r) { if (r >= H) return 0; return 40 / (Math.PI * H2 * H2 * H2) * (H - r); }

const REST_DENSITY = 300; // target density (kg/m²-ish, in particle units) const MASS = 1; // particle mass function computeDensities(particles) { for (const pi of particles) { let density = 0; for (const pj of particles) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r2 = dx*dx + dy*dy; density += MASS * Wpoly6(r2); } pi.density = Math.max(density, REST_DENSITY * 0.01); // avoid zero } }

const STIFFNESS = 200; // pressure = stiffness * (density - rest_density) function computePressure(p) { p.pressure = STIFFNESS * (p.density - REST_DENSITY); } function pressureForce(pi, pj) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r = Math.sqrt(dx*dx + dy*dy); if (r === 0 || r >= H) return { x: 0, y: 0 }; // Symmetric pressure force (Newton's 3rd law) const pressure = (pi.pressure + pj.pressure) / 2; const grad = gradWspiky(dx, dy, r); return { x: -MASS * pressure / pj.density * grad.x, y: -MASS * pressure / pj.density * grad.y }; }

const VISCOSITY = 200; function viscosityForce(pi, pj) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r = Math.sqrt(dx*dx + dy*dy); if (r >= H) return { x: 0, y: 0 }; const lap = lapWviscosity(r); const dvx = pj.vx - pi.vx; const dvy = pj.vy - pi.vy; return { x: VISCOSITY * MASS / pj.density * dvx * lap, y: VISCOSITY * MASS / pj.density * dvy * lap }; } // Combine all forces function computeForces(particles, gravity = 500) { for (const pi of particles) { pi.fx = 0; pi.fy = 0; computePressure(pi); } for (let i = 0; i < particles.length; i++) { const pi = particles[i]; for (let j = i + 1; j < particles.length; j++) { const pj = particles[j]; const pf = pressureForce(pi, pj); const vf = viscosityForce(pi, pj); pi.fx -= pf.x + vf.x; pi.fy -= pf.y + vf.y; pj.fx += pf.x + vf.x; pj.fy += pf.y + vf.y; } pi.fy += gravity * pi.density; // gravity proportional to density } }

const DT = 0.003, DAMPING = -0.5; const W = 800, HEIGHT = 600; // canvas size function integrate(particles) { for (const p of particles) { // Semi-implicit (symplectic) Euler p.vx += DT * p.fx / p.density; p.vy += DT * p.fy / p.density; p.x += DT * p.vx; p.y += DT * p.vy; // Boundary collisions if (p.x < H) { p.x = H; p.vx *= DAMPING; } if (p.x > W - H) { p.x = W - H; p.vx *= DAMPING; } if (p.y < H) { p.y = H; p.vy *= DAMPING; } if (p.y > HEIGHT - H) { p.y = HEIGHT - H; p.vy *= DAMPING; } } } // Render function render(ctx, particles) { ctx.clearRect(0, 0, W, HEIGHT); ctx.fillStyle = 'rgba(0,10,30,1)'; ctx.fillRect(0, 0, W, HEIGHT); for (const p of particles) { const speed = Math.sqrt(p.vx*p.vx + p.vy*p.vy); const t = Math.min(speed / 300, 1); ctx.fillStyle = `hsl(${200 - t*140}, 80%, ${40 + t*30}%)`; ctx.beginPath(); ctx.arc(p.x, p.y, H * 0.3, 0, Math.PI*2); ctx.fill(); } }

class SpatialHash { constructor(cellSize) { this.cs = cellSize; this.table = new Map(); } key(x, y) { return `${Math.floor(x/this.cs)},${Math.floor(y/this.cs)}`; } clear() { this.table.clear(); } insert(p) { const k = this.key(p.x, p.y); if (!this.table.has(k)) this.table.set(k, []); this.table.get(k).push(p); } neighbours(p, radius) { const result = []; const cells = Math.ceil(radius / this.cs); const cx = Math.floor(p.x / this.cs), cy = Math.floor(p.y / this.cs); for (let dx = -cells; dx <= cells; dx++) for (let dy = -cells; dy <= cells; dy++) { const bucket = this.table.get(`${cx+dx},${cy+dy}`); if (bucket) result.push(...bucket); } return result; } } // Use in computeForces: replace inner loop with // const neighbours = spatialHash.neighbours(pi, H); // for (const pj of neighbours) { ... } // Time complexity: O(N * k) where k = avg neighbours (~20–50)