Introduction

In 1977, two independent teams — Gingold & Monaghan and Lucy — needed to simulate astrophysical gas clouds without a fixed computational grid. Grids struggle when matter tears apart or collapses; the teams' solution was radical: make the fluid itself carry the equations. The result, Smoothed-Particle Hydrodynamics (SPH), is still in active use half a century later.

The core idea is deceptively simple. Instead of carving space into cells, you scatter particles through the fluid. Each particle carries its own mass, velocity, density and pressure. To evaluate any quantity at a point in space, you sum the contributions of nearby particles, each weighted by a smoothing kernel — a bell-shaped function W(r,h)W(r, h) that falls to zero beyond a radius hh called the smoothing length.

Because the particles move with the fluid, SPH is Lagrangian: it tracks fluid elements rather than fixed points in space. There is no mesh to tangle, no cells to refill — just particles following their own trajectories while continuously averaging over their neighbors.

This article explains how that averaging works, what it costs, and why SPH shows up in everything from numerical simulations of colliding stars to the water in your favorite blockbuster film.

Try It: Watch a Splash

The demo below runs a miniature SPH simulation entirely in your browser. A blob of particles drops under gravity, hits the floor, and settles. Each particle pushes its neighbors through a pressure force proportional to the density gradient computed by the smoothing kernel.

<!-- SPH splash demo — HTML structure -->
<div class="controls">
  <label>Gravity: <input type="range" id="gravity" min="100" max="800" value="400" step="10"></label>
  <label>Damping: <input type="range" id="damping" min="0" max="99" value="20" step="1"></label>
  <button id="btn-drop" type="button">Drop particles</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>
<canvas id="c" width="480" height="340"></canvas>
<div id="status" class="status">Press "Drop particles" to start the simulation.</div>

/* SPH demo styles — canvas + controls */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; background: #f0f4f8; }
.controls { display: flex; gap: .6rem; flex-wrap: wrap; align-items: center;
            padding: .5rem .4rem; background: #e2eaf2; border-radius: 8px; margin-bottom: .5rem; }
label { font-size: .82rem; color: #334; display: flex; align-items: center; gap: .3rem; }
input[type=range] { width: 80px; cursor: pointer; }
canvas { display: block; border-radius: 8px; border: 1px solid #c4d0dc; background: #dce8f5; width: 100%; }
.status { font-size: .88rem; font-weight: 600; color: #2a4d6e; margin-top: .4rem; min-height: 1.3em; }
button { font: 600 13px system-ui; padding: .38rem .8rem; border: 1px solid #2a4d6e;
         background: #2a4d6e; color: #fff; border-radius: 6px; cursor: pointer; }
button.ghost { background: #fff; color: #2a4d6e; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* SPH (Smoothed-Particle Hydrodynamics): fluid is represented as particles. */
/* Each particle samples density and pressure from its neighbors via a kernel. */
const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

/* Simulation parameters */
const NUM = 80;
const RADIUS = 28;       // smoothing radius in pixels
const REST_DENSITY = 8;  // target rest density
const STIFFNESS = 800;   // pressure stiffness constant
const VISCOSITY = 0.3;   // viscosity coefficient
const PARTICLE_R = 5;    // visual radius

let gravity = 400;
let damping = 0.20;
let animId = null;
let running = false;
let particles = [];

/* Poly6 smoothing kernel for density estimation */
function poly6W(r2, h) {
  const h2 = h * h;
  if (r2 >= h2) return 0;
  const d = h2 - r2;
  return (315 / (64 * Math.PI * Math.pow(h, 9))) * d * d * d;
}

/* Spiky gradient kernel for pressure force */
function spikyGrad(r, h) {
  if (r <= 0 || r >= h) return 0;
  const d = h - r;
  return -(45 / (Math.PI * Math.pow(h, 6))) * d * d;
}

/* Viscosity Laplacian kernel */
function viscLaplacian(r, h) {
  if (r >= h) return 0;
  return (45 / (Math.PI * Math.pow(h, 6))) * (h - r);
}

function initParticles() {
  particles = [];
  /* Drop particles in a rectangular blob near top-center */
  const cols = 8, rows = Math.ceil(NUM / cols);
  const sx = W * 0.35, sy = H * 0.05;
  const spacing = PARTICLE_R * 2.2;
  for (let i = 0; i < NUM; i++) {
    const col = i % cols, row = Math.floor(i / cols);
    particles.push({
      x: sx + col * spacing + (Math.random() - 0.5) * 2,
      y: sy + row * spacing + (Math.random() - 0.5) * 2,
      vx: 0, vy: 0,
      density: REST_DENSITY,
      pressure: 0
    });
  }
}

function computeDensityPressure() {
  /* Estimate density at each particle by summing kernel contributions */
  for (const pi of particles) {
    pi.density = 0;
    for (const pj of particles) {
      const dx = pj.x - pi.x, dy = pj.y - pi.y;
      const r2 = dx * dx + dy * dy;
      pi.density += poly6W(r2, RADIUS);
    }
    pi.density = Math.max(pi.density, 0.001);
    pi.pressure = STIFFNESS * (pi.density - REST_DENSITY);
  }
}

function computeForces(dt) {
  /* Compute pressure + viscosity forces and integrate positions */
  for (const pi of particles) {
    let fx = 0, fy = gravity;
    for (const pj of particles) {
      if (pi === pj) continue;
      const dx = pj.x - pi.x, dy = pj.y - pi.y;
      const r = Math.sqrt(dx * dx + dy * dy);
      if (r < 0.001 || r >= RADIUS) continue;

      /* Pressure gradient force (symmetric form) */
      const gradW = spikyGrad(r, RADIUS);
      const pterm = -(pi.pressure + pj.pressure) / (2 * pj.density) * gradW / r;
      fx += pterm * dx;
      fy += pterm * dy;

      /* Viscosity diffusion force */
      const lapW = viscLaplacian(r, RADIUS);
      fx += VISCOSITY * (pj.vx - pi.vx) / pj.density * lapW;
      fy += VISCOSITY * (pj.vy - pi.vy) / pj.density * lapW;
    }
    /* Semi-implicit Euler integration */
    pi.vx += fx * dt / pi.density;
    pi.vy += fy * dt / pi.density;

    /* Boundary collision with energy loss (damping) */
    const margin = PARTICLE_R;
    if (pi.x < margin) { pi.x = margin; pi.vx = Math.abs(pi.vx) * (1 - damping); }
    if (pi.x > W - margin) { pi.x = W - margin; pi.vx = -Math.abs(pi.vx) * (1 - damping); }
    if (pi.y < margin) { pi.y = margin; pi.vy = Math.abs(pi.vy) * (1 - damping); }
    if (pi.y > H - margin) { pi.y = H - margin; pi.vy = -Math.abs(pi.vy) * (1 - damping); }

    pi.x += pi.vx * dt;
    pi.y += pi.vy * dt;
  }
}

function draw() {
  ctx.clearRect(0, 0, W, H);
  /* Color by speed — blue=slow, red=fast */
  for (const p of particles) {
    const speed = Math.sqrt(p.vx * p.vx + p.vy * p.vy);
    const t = Math.min(speed / 200, 1);
    const r = Math.round(40 + t * 180), b = Math.round(200 - t * 140);
    ctx.beginPath();
    ctx.arc(p.x, p.y, PARTICLE_R, 0, Math.PI * 2);
    ctx.fillStyle = `rgba(${r},120,${b},0.85)`;
    ctx.fill();
  }
}

let last = null;
function loop(ts) {
  if (!running) return;
  if (last === null) last = ts;
  const dt = Math.min((ts - last) / 1000, 0.02);
  last = ts;
  computeDensityPressure();
  computeForces(dt);
  draw();
  animId = requestAnimationFrame(loop);
}

function startDrop() {
  if (running) { running = false; cancelAnimationFrame(animId); last = null; }
  initParticles();
  draw();
  running = true;
  last = null;
  animId = requestAnimationFrame(loop);
  document.getElementById('status').textContent = i18nFmt('Simulating %n% SPH particles — watch them splash and settle.', { n: NUM });
}

function reset() {
  running = false;
  cancelAnimationFrame(animId);
  last = null;
  initParticles();
  draw();
  document.getElementById('status').textContent = 'Press "Drop particles" to start the simulation.';
}

document.getElementById('gravity').addEventListener('input', function() { gravity = +this.value; });
document.getElementById('damping').addEventListener('input', function() { damping = +this.value / 100; });
document.getElementById('btn-drop').addEventListener('click', startDrop);
document.getElementById('btn-reset').addEventListener('click', reset);

initParticles();
draw();
document.getElementById('status').textContent = 'Press "Drop particles" to start the simulation.';
Result
Console

Try adjusting gravity and damping with the sliders. Notice how increasing the smoothing radius makes the fluid look more viscous — particles "feel" more neighbors and average out sudden pressure spikes. Decrease it too much and the simulation becomes noisy: each particle sees too few neighbors and pressure estimates jump wildly. The sweet spot is roughly 30–50 neighbors per particle, a rule of thumb SPH practitioners have used since the 1980s.

The Real Complexity

SPH is elegant, but naive implementations hide a cost that scales badly.

One time step breaks into three phases:

  1. Neighbor search — for each particle, find every other particle within radius hh.
  2. Density estimation — sum ρi=∑jmjW(∣ri−rj∣,h)\rho_i = \sum_j m_j W(|\mathbf{r}_i - \mathbf{r}_j|, h) over neighbors.
  3. Force evaluation — apply the SPH form of the Navier-Stokes equations: pressure gradient, viscosity and any external forces.

With NN particles, the naive pairwise search costs O(N2)O(N^2) per step — fine for a few hundred particles, disastrous at a million.

The fix: spatial hashing. Divide space into cells of size hh. Each particle hashes to its cell; to find neighbors you only look at the 3d3^d adjacent cells (32=93^2 = 9 in 2-D, 33=273^3 = 27 in 3-D). If particles are roughly uniformly spread, the neighbor count per particle is bounded by a constant, dropping the per-step cost to O(N)O(N) for force evaluation. The hash table rebuild costs O(N)O(N) with a counting sort, giving an overall per-step complexity of O(Nlog⁡N)O(N \log N) in practice (dominated by sort on GPU) or O(N)O(N) in ideal uniform distributions.

Accuracy scales with both NN (more particles) and hh (smaller kernel). Reducing hh improves spatial resolution but forces smaller time steps (the CFL condition: Δt≤C⋅h/vmax⁡\Delta t \leq C \cdot h / v_{\max}), so the total cost to reach time TT scales roughly as O(N4/3)O(N^{4/3}) in three dimensions for fixed accuracy — an uncomfortable truth for high-resolution ocean simulations.

SPH sits in a crowded zoo of meshfree methods alongside the Material Point Method (MPM) and Moving Least Squares approaches. Each trades accuracy, stability and cost differently, but the neighbor-search bottleneck is universal: no matter how you smooth, you must still find who is close to whom.

Where It Matters

Because SPH needs no grid, it excels wherever geometry changes violently:

  • Astrophysics: the method's birthplace. SPH simulates galaxy mergers, proto-planetary disk formation and supernova shockwaves. The GADGET and AREPO codes run SPH at cosmological scales.
  • Tsunami and dam-break modeling: water that tears across irregular terrain is exactly the free-surface, large-deformation regime where grids fail and SPH thrives.
  • Metal forming and impact mechanics: solid materials undergoing fracture and plastic flow are treated as nearly-incompressible fluids in SPH, modeling bullet impacts, explosions and crash tests.
  • Visual effects: DreamWorks, ILM and Sony Pictures Imageworks have used SPH-based solvers (often coupled with FLIP) to animate water, lava and crowds of tumbling objects in films and games.
  • GPU acceleration: SPH maps cleanly onto massively parallel hardware. NVIDIA's PhysX and Bullet Physics both include GPU SPH for real-time fluid in games.

The same kernel-interpolation idea that makes SPH work for fluids also appears in kernel density estimation in statistics and in radial basis function interpolation — a reminder that the mathematical core of SPH is more general than any single application domain. It connects naturally to the broader landscape of numerical simulations and optimization problems that drive scientific computing.

Conclusion

SPH's bet is simple: if the fluid is the solver, there is no mesh to corrupt, no cell to flood, no topology to repair. Particles carry mass, momentum and the smoothing kernel does the rest — integrating the Navier-Stokes equations one neighbor-weighted sum at a time.

That bet has paid off across fifty years and a dozen disciplines. From the galaxies Gingold and Monaghan wanted to collide in 1977, to the oceans of the latest blockbuster, SPH's meshfree Lagrangian approach keeps proving itself wherever matter refuses to stay still.

The price is the neighbor search: as NN grows, keeping that search efficient is what separates a toy simulation from a production fluid engine. Master the spatial hash, tune the smoothing length, and you hold one of the most versatile simulation tools in computational science.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/smoothed-particle-hydrodynamics/Content licensed under CC BY-NC 4.0.