How to Debug a SPH Fluid Simulation: My Step-by-Step Playbook
Marcel Kießlich·
My first SPH simulation didn't explode dramatically. It just... showed nothing. A blank canvas, a console full of NaN, and no clue where the numbers had gone wrong. I'd stare at the code, convince myself the math was right, run it again — NaN. For weeks, every run was either a blank screen or particles shooting off to infinity, and I couldn't tell which part of the pipeline was broken.
When I finally killed the NaN bugs, the reward was a different kind of wrong: a swirling, gassy mess that looked more like smoke than water. Turns out I was using the basic linear equation of state and had no viscosity at all. The particles were technically alive, just behaving like a gas in a box.
If you've written your own Smoothed Particle Hydrodynamics solver, you know these exact stages. SPH is deceptively easy to get almost right. The math fits on a napkin. The hard part is that a single wrong sign, a slightly-too-large time step, or a neighbor list quietly dropping particles all produce the same useless result: garbage on screen.
So instead of another wall of theory, here's the actual order I work in when I'm building or debugging a 2D SPH fluid from scratch. This playbook has saved me more time than any paper. Below it you'll find the deeper explanations for each step.
The playbook: how I debug SPH, in order
1. Always build it in 2D first. The math and the bugs are identical to 3D, but in 2D everything is faster to run, trivial to visualize, and you can actually watch what each particle is doing. Port to 3D only once 2D is rock solid. Debugging a blowup in a 3D point cloud is masochism; debugging it in a flat 2D box is a Tuesday.
2. Rule out the neighbor search before anything else. A broken neighbor search makes everything explode, because density, pressure, and forces are all neighbor sums — so a bug here masquerades as a dozen other bugs. Temporarily rip out your spatial grid and replace it with a brute-force "check every pair" loop. It's slow, but it's almost impossible to get wrong. If the brute-force version is stable and the grid version explodes, you've found your bug in one move. Only optimize the neighbor search back in once the physics is proven.
3. Turn gravity off and spawn the particles in the middle of empty space. This is my favorite isolation test. With no gravity and no walls, a blob of particles sitting in the middle of the domain should gently relax to its rest spacing and then basically sit there. If it quietly settles, your core pressure/density loop is healthy. If it explodes outward or implodes into a clump, you've caught a force bug with every external influence stripped away. Add gravity and boundaries back only after this passes.
4. Log everything and make sure the numbers are sane. Print density and pressure for a few particles every step. Interior density should hover right around your rest density ρ0; pressure should be small and positive near rest, larger under compression. If an interior particle reports a density of 3 when rest is 1000, or a pressure of 1012, stop — you've found the problem before it even reaches the force step. Numbers lie less than pixels do.
5. Clamp pressure to zero.p = max(0, ...). Near a free surface, particles have neighbors on only one side, so density reads low and your equation of state computes negative pressure, which physically means "suck inward" — and that's what makes a fluid curdle into clumps. Refusing negative pressure is the single highest-impact one-liner in this whole article.
6. Don't use the textbook gas law — use the Tait equation. The simple linear p=k(ρ−ρ0) is springy and makes your liquid look weirdly gassy and bouncy. Switch to Tait (details and code below). The 7th-power term barely resists tiny compressions but slams back hard against big ones, which is exactly what reads as an incompressible liquid.
7. If your fluid looks like a gas and swirls around, add viscosity — I'd use XSPH. A noisy, swirling, never-settling fluid usually has too little damping of particle-scale noise. Raw artificial viscosity works but dumps energy out of the dynamics. XSPH velocity smoothing calms the chaos by nudging each particle toward its neighbors' average velocity for advection only, without poisoning the momentum equation. Formula below.
8. Use the right kernels — and know that 2D kernels are not 3D kernels. The normalization constants in Müller's poly6, spiky, and viscosity kernels are derived for 3D. Paste those straight into a 2D sim and your densities and forces are silently mis-scaled. The 2D constants (popularized by Solenthaler et al.'s SPH Based Shallow Water Simulation) are different. Exact values below — get these right and a lot of "mysterious" misbehavior just disappears.
That's the loop. The rest of this post is the why and the how behind each step.
A 30-second mental model
You can't debug SPH without a gut feel for what it does, so here's the whole method in a paragraph.
Every particle carries mass and properties. To compute anything — density, pressure force, viscosity — you sum contributions from neighbors inside the support radius h, weighted by a smoothing kernel W(r,h). Density is a sum of neighbor masses times the kernel. Pressure comes from density via an equation of state. Pressure force comes from the gradient of the kernel. That's it. Almost every bug is a corruption of one of those four things, so keep the chain in your head: neighbors → density → pressure → force. When something blows up, the question is always which link broke — which is exactly what steps 2 through 4 of the playbook are designed to localize.
"My particles explode" — the classic blowup
The most common first failure. The sim looks fine, then energy appears from nowhere and it detonates. Four usual suspects, roughly in order:
Your time step is too big
Weakly compressible SPH uses a high artificial sound speed to fake incompressibility, and that fast wave speed bounds your step via the CFL condition. A reliable starting set of constraints:
cs≈10⋅vmax
Δtcfl=0.4csh
Δtforce=0.25∣amax∣h
Δtvisc=0.4νh2
Δt=min(Δtcfl,Δtforce,Δtvisc)
If you halve Δt and the explosion delays or vanishes, this was it. It often is.
Your equation of state is too stiff
Crank stiffness up to enforce incompressibility and the system gets so stiff any reasonable step blows up. Tune stiffness and Δt together, never separately.
Your pressure force isn't symmetric
If particle i pushes on j harder than j pushes back, you inject momentum every step and it accumulates into a delayed eruption. Use a momentum-conserving form:
fi+=−mimj(ρi2pi+ρ∇Wspiky(rij,h)
Your kernel gradient has the wrong sign
Flipping it turns repulsion into attraction. The tell is an immediate, violent explosion rather than a slow drift. Sanity check: two slightly-compressed particles must push apart.
"My particles clump into clusters" — tensile instability and pairing
Instead of smooth fluid, particles gather into tight clumps with empty gaps, like the water curdled. Two causes — and they have different names in the literature.
Negative pressure (tensile instability), exactly as in playbook step 5. Near a free surface, density drops below rest, the equation of state returns negative pressure, and particles get sucked together. Clamp pressure to zero and most of this vanishes.
The wrong kernel for the pressure force (pairing instability). This one is more subtle. Use poly6 for density but the spiky kernel for the pressure gradient. The poly6 gradient flattens to zero as particles get very close (∇Wpoly6∝−6r(h2−r2)2→0 as r→0), so two near-coincident particles barely repel and happily merge — they "pair up." The spiky gradient stays strong near the center and keeps them apart. Swap these two and you'll lose an afternoon to a bug that isn't in your force code at all.
"Everything turns into NaN"
One frame the numbers are real, the next they're NaN, and from then on it spreads through every sum.
The dominant cause is two particles at exactly the same position: distance r is zero, you normalize the direction by dividing by r, and you've divided by zero.
float r = length(r_ij);
if (r < 1e-6f) continue; // skip self / coincident particles
vec2 dir = r_ij / r;
The other cause is subtler but just as deadly: density hitting zero. Look at the pressure force equation:
fi+=−mimj(ρi2pi+ρ∇Wspiky(rij,h)
See the ρi2 and ρj2 in the denominators? If any particle's density is zero, you divide by zero and everything is NaN from that frame onward. This happens when a particle has no neighbors within the support radius — its density sum comes out to zero because there's nothing contributing to it.
The fix: always include each particle's own self-contribution to density. Every particle contributes to its own density via W(0,h), which is simply the kernel evaluated at zero distance:
ρi=miW(0,h)+∑j=imjW(rij,h)
The first term is the self-density. Even if a particle is completely alone in empty space with zero neighbors, its density will still be mi⋅W(0,h) — a positive, finite number. Division by zero is impossible.
// Initialize density with self-contribution BEFORE the neighbor loopfor (int i = 0; i < n; i++) {
rho[i] = mass * W_poly6(0.0f, h); // never zero
}
// Then accumulate neighbor contributionsfor (each neighbor j of i) {
float r = length(pos[i] - pos[j]);
rho[i] += mass * W_poly6(r, h);
}
Forgetting this self-contribution is one of the most common beginner mistakes — it looks harmless because it "just" adds a constant to every density, but without it you're one isolated particle away from a NaN cascade. Add an assertion that fires the instant any value goes non-finite and prints the particle index, because by the time you see NaNs on screen you're far from the source.
Equation of state: use Tait, not the textbook gas law
This is playbook step 6, in detail. Most tutorials show the linear gas law:
p=k(ρ−ρ0)
It's linear, so under large compression it doesn't push back hard enough, and the fluid reads as bouncy and gas-like. The fix is the Tait equation (weakly compressible SPH, Becker & Teschner 2007):
p=B[(ρ0ρ)γ−1]
where:
γ=7,B=γρ0cs2
and cs is the artificial speed of sound, typically cs≈10⋅vmax.
The 7th power barely resists small compressions but resists large ones extremely hard, which is precisely what makes the fluid behave like an incompressible liquid instead of a springy gas. You can still clamp the result to zero for calmer surfaces:
float p = fmax(0.0f, B * (powf(rho / rho0, 7.0f) - 1.0f));
If your fluid swirls like a gas: add viscosity (use XSPH)
Playbook step 7. A fluid that never settles, swirls endlessly, and looks more like smoke than water usually has too little damping of particle-scale noise. You can add classic artificial viscosity, but it dumps real energy out of the dynamics and can over-damp into jelly. I prefer XSPH velocity smoothing (Monaghan), which nudges each particle's velocity toward the neighborhood average and uses that smoothed velocity only to move the particle:
vismooth=vi+ε∑jρjmj(vvi)W(rij,h)
xi+=Δt⋅vismooth
with ε≈0.1–0.5.
This kills the swirl and the noise without injecting viscosity into the momentum equation, so the fluid stays lively but coherent. If you genuinely want a thick, syrupy fluid, then reach for a real viscosity term — but for "stop looking like a gas," XSPH is the cleaner tool.
Get your 2D kernels right (they are not the 3D ones)
Playbook step 8, the one people skip and regret. Müller's poly6, spiky, and viscosity kernels are normalized for 3D. Drop those constants into a 2D simulation and your density, viscosity, and external-force scaling are all silently wrong. The 2D-renormalized versions — the set commonly cited from Solenthaler et al., SPH Based Shallow Water Simulation — are what you want.
Density (poly6) — 2D
Wpoly6(r,h)=πh84(h2−r2)3,0≤r≤h
Pressure (spiky gradient) — 2D
∇Wspiky(r,h)=−πh510(h−r)2∣r∣r
Note: Strictly differentiating and normalizing the spiky kernel in 2D gives −30/(πh5), not −10/(πh5). Both values appear in the literature and both work — the difference gets absorbed into your pressure stiffness constant. Just make sure your code and your mental model use the same one.
Viscosity (Laplacian) — 2D
∇2Wvisc(r,h)=πh540(h−r)
For contrast, the 3D constants are completely different:
Kernel
3D constant
2D constant
Poly6
64πh9315
πh84
Spiky gradient
−πh645
−πh
Viscosity Laplacian
πh645
πh
Same kernel shapes, different normalization — because you're integrating over a disk in 2D and a ball in 3D.
One honest note: pressure and surface tension don't strictly need a perfectly normalized constant, because the factor cancels out in those computations and effectively gets absorbed into your stiffness tuning. Density, viscosity, and external forces do care. So if you must rush, get the poly6 (density) constant exactly right first.
"Particles leak through the walls" — boundary problems
Fast particles tunnel through solid walls in a single step before any penalty force can act — classic tunneling. Fixes: smaller time steps, a collision response that clamps position and reflects the velocity component normal to the wall (reflecting position but not velocity leaves the particle jittering into the wall forever), or — most robustly — boundary particles: a line of stationary particles lining the walls that participate in the density and pressure sums exactly like fluid. Real particles feel the wall's pressure before reaching it, which also fixes the density underestimate near boundaries.
"The fluid behaves like jelly" — too much damping
The opposite of exploding: sluggish, refuses to splash, surface goes flat instantly. Almost always too much artificial viscosity. Turn it down until you see noise, then back off slightly. And check you're not double-damping — artificial viscosity plus a physical viscosity term plus a velocity decay on integration all stack into mush. This is exactly why I prefer XSPH (above) for visual smoothness without energy loss.
Frequently asked questions
Why do my SPH particles explode after a few frames?
Usually a time step that violates the CFL condition, an equation of state that's too stiff, or a non-symmetric pressure force injecting momentum. Halve your time step first; if the explosion delays or vanishes, that was it. If it doesn't, rule out the neighbor search with a brute-force loop.
Why are my SPH particles clumping together?
Tensile instability — usually negative pressure near the surface (clamp with p = max(0, ...)) or using the poly6 gradient for the pressure force instead of the spiky kernel.
What equation of state should I use for SPH water?
The Tait equation: p=B[(ρ/ρ0)7−1] with B=ρ0cs2/7. The linear gas law makes liquids look springy and gas-like; the 7th-power Tait term reads as a proper incompressible liquid.
My SPH fluid swirls like a gas and never settles — how do I fix it?
Add velocity smoothing. XSPH (vsmooth=v+ε∑jρjmj(vj−vi)W) calms the swirl without dumping energy out of the dynamics like raw artificial viscosity does.
Are 2D and 3D SPH kernels the same?
No. The kernel shapes match but the normalization constants differ because you integrate over a disk versus a ball. In 2D, poly6 uses 4/(πh8), the spiky gradient uses −10/(πh5), and the viscosity Laplacian uses 40/(πh5). Using 3D constants in a 2D sim silently mis-scales everything.
Why is my SPH simulation producing NaN values?
Two coincident particles cause a divide-by-zero when normalizing the direction vector. Skip pairs closer than a small epsilon, and include each particle's self-contribution so density never underflows to zero.
The one habit that ties it all together
If you take away one thing: build your visual and logging debug tools before you build features. Color-by-density, color-by-neighbor-count, a force-vector overlay, a single-step button, a NaN tripwire, and console logs of density and pressure. Combined with the playbook — 2D first, brute-force the neighbor search, gravity off and spawn in the middle, sanity-check the numbers — almost every bug in this article stops being mysterious and starts being something you can see.
SPH isn't fragile because the method is bad. It's fragile because everything is coupled, so a small error anywhere surfaces as the same generic catastrophe. Once you can localize which link in the chain broke, debugging it becomes almost mechanical. Almost.