In a previous post, we simulated sloshing in a rectangular tank using CLSVOF on a fixed grid. It worked, but only barely interactively, at coarse resolution, and only for a box. The grid made arbitrary container shapes impractical, and the visual quality at interactive resolutions was limited. This time, we wanted to see if a meshless approach could do better: higher visual quality, arbitrary containers, and genuinely interactive frame rates. No grid — just particles.

The container is a Stanford bunny — non-convex, thin-featured, immediately recognizable. Fill it with ~20,000 SPH particles, let the user grab and shake it, and render the result as water in glass. The simulation runs at about 20 FPS on an Apple M3 laptop. Not silky smooth, but responsive enough to play with. This post covers the physics, the rendering, and where we ran into the fundamental limits of purely particle-based fluid simulation.

The constraint was purely Lagrangian: no background grid, no hybrid methods. Just particles tracked through space. This is the SPH (Smoothed Particle Hydrodynamics) family of methods, where the fluid is discretized into particles that carry mass, velocity, and pressure, and interact through kernel functions.

How SPH works

In SPH, every fluid quantity — density, pressure, velocity — is interpolated from nearby particles using a smooth kernel function W. A particle’s density, for example, is the weighted sum over its neighbors:

\rho_i = \sum_j m_j \, W(\mathbf{x}_i – \mathbf{x}_j, h)

where h is the smoothing radius that defines how far each particle’s influence reaches. Gradients of these fields are computed analytically from the kernel gradient \nabla W, which is what makes pressure forces, viscosity, and incompressibility correction possible without any mesh or grid.

The kernel functions matter. We use two: the poly6 kernel for density estimation (smooth, cheap to evaluate) and the spiky kernel for pressure gradients (its sharp peak near r = 0 prevents particles from collapsing onto each other).

Neighbor search is the practical bottleneck. Every particle needs to know who’s within distance h. We use a spatial hash — a uniform grid of cell size h, where each particle hashes into a cell. Querying a particle’s neighborhood means checking its cell and the 26 surrounding cells. With a counting-sort layout (CSR), the neighbor lists are contiguous in memory, which matters when you’re doing this 20,000 times per frame.

Choosing a solver: PBF, then DFSPH

We started with Position Based Fluids (PBF, Macklin & Muller 2013). PBF enforces incompressibility by iteratively projecting particle positions to satisfy a density constraint — similar to how position-based dynamics handles cloth or rigid bodies. It’s simple, stable under large deformations, and converges in a few iterations.

PBF worked well enough to get the simulation running. But at the iteration counts we could afford in real time (~3-5 per substep), particles never quite settled. They retained a residual oscillation of about 0.15 m/s — visibly vibrating even when the container was stationary. We tried Gauss-Seidel ordering, higher viscosity, divergence correction, and temporal damping. Some helped, but none eliminated the jitter within our performance budget. PBF can converge given enough iterations, but “enough” was more than we could spend.

Divergence-free SPH (DFSPH, Bender & Koschier 2015) converges faster for the same iteration budget. It corrects velocities instead of positions, using two iterative loops:

1. Constant density correction — adjusts velocities so that the predicted density matches the rest density.
2. Divergence-free correction — removes velocity divergence, preventing local compression before it happens.

Both loops use precomputed alpha factors — essentially the inverse of how strongly each particle’s neighbors respond to velocity changes. The alpha factor for particle i is:

\alpha_i = \frac{\rho_i}{\left|\sum_j \nabla W_{ij}\right|^2 + \sum_j \left|\nabla W_{ij}\right|^2}

This makes each correction step well-scaled regardless of local particle configuration. The alpha factor is folded into a per-particle stiffness that encodes the density error,

\kappa_i = \frac{\alpha_i \, (\rho_i^* – \rho_0)}{\Delta t^2}

and the velocity correction per iteration is:

\Delta \mathbf{v}_i = -\Delta t \sum_j \left(\kappa_i + \kappa_j\right) \nabla W_{ij}

where \rho_i^* is the predicted density. Convergence is fast — typically 3-5 iterations per loop — and the result is smooth, stable fluid that actually comes to rest.

The rest density trap

One thing that tripped us up early: the rest density \rho_0. In textbook SPH, this is typically 1000 kg/m³ for water. But our solver uses raw kernel sums without explicit particle masses, so the numerical density at equilibrium spacing depends entirely on the kernel shape and particle spacing. Setting \rho_0 = 1000 when the actual kernel-evaluated density was around 7500 caused the solver to violently push all particles apart.

The fix: auto-compute \rho_0 by sampling the kernel density at a few interior particles after seeding. This makes the solver robust to any particle radius and smoothing length without manual tuning. A small thing, but it would have been completely baffling without understanding that the density constraint operates on kernel sums, not physical densities.

The container: signed distance fields

The bunny mesh has ~28,000 triangles. Testing every particle against every triangle each frame is out of the question. Instead, we voxelize the mesh into a signed distance field (SDF) — a 128^3 grid where each voxel stores the distance to the nearest surface, negative inside and positive outside.

Building the SDF uses a BVH (bounding volume hierarchy) for fast closest-point queries. Inside/outside classification uses a majority-vote ray test: cast three rays in different directions, count intersections, and take the majority. This is more robust than a single ray for non-watertight or self-intersecting meshes.

At runtime, any point’s distance to the container wall is a single trilinear interpolation — O(1). The gradient of the SDF (via central differences) gives the wall normal for collision response. Particles that drift outside are reprojected onto the surface with velocity reflection and friction.

A key architectural choice: the simulation runs in the container’s local coordinate frame. When the user rotates the container, we rotate the gravity vector, not the particles. This means the SDF never needs recomputation.

Getting a water surface: marching cubes and Zhu-Bridson

Particles alone don’t give you a water surface. To render the fluid as a continuous sheet of water, we need to extract an iso-surface — and that’s where Marching Cubes comes in.

The classic approach is to evaluate the SPH density field on a regular grid and extract the surface where density equals the rest density. We implemented this, and it mostly worked — except near container walls.

The wall problem

SPH density estimation breaks down at boundaries. A particle near a wall has neighbors only on one side, so its kernel sum drops to roughly 50% of what an interior particle sees. The density field sags near every wall, producing curved, inward-pulled surfaces where the water should meet the glass cleanly.

We tried several fixes: ghost particles mirrored outside the container, density clamping, correction factors. None worked well. The fundamental issue is that density-based surface reconstruction can’t distinguish “near a wall” from “near the free surface.”

The solution was Zhu-Bridson particle level sets (Zhu & Bridson 2005). Instead of evaluating density, for each grid point you:

1. Find nearby particles
2. Compute their kernel-weighted average position (centroid)
3. Set \phi = r_p – |\mathbf{x} – \mathbf{x}_{\text{centroid}}|

where r_p is the particle radius. Positive inside the fluid, negative outside. This is a geometric measure — distance to the local particle cloud center — so it doesn’t suffer from the one-sided kernel problem at walls. Combined with clamping against the container SDF (\phi = \min(\phi_{\text{fluid}}, \phi_{\text{container}})), the water surface meets the glass cleanly.

Mesh post-processing

The raw Marching Cubes output is jagged and needs post-processing. We clamp vertices near the container wall onto the SDF surface (creating a clean waterline), then apply Laplacian smoothing and a midpoint subdivision pass. We found it more effective to smooth the coarse mesh before subdividing — after subdivision the edges are short enough that the smoothing has little effect. The Laplacian we use is the simple combinatorial one (uniform neighbor weights), which isn’t geometrically correct but produces acceptable results for visual purposes. A cotangent Laplacian would be more principled but requires a Delaunay triangulation, which MC output doesn’t guarantee.

Wall vertices get their normals from the container SDF gradient rather than from the triangle geometry, which gives a crisp join between water and glass. A temporal blend (mixing 80% of the current field with 20% of the previous frame) suppresses frame-to-frame jitter in the surface, and lets us skip the Zhu-Bridson evaluation on alternate frames — saving ~8ms per frame with almost no visible difference.

Rendering

The rendering uses standard real-time water techniques: the scene is rendered to an off-screen buffer for refraction, a caustic map is generated by displacing water mesh vertices along refracted light rays onto the container floor (Wyman 2005), and the water surface is composited with Fresnel-driven reflection/refraction and Beer’s law absorption. The glass container samples the same HDR environment cubemap for reflections. Particles classified as spray (fewer than 12 neighbors) are rendered separately as small shaded spheres.

Optimization: from 8 to 20 FPS

With ~20,000 particles, the initial implementation ran at about 8 FPS — too slow to feel interactive. After optimization it’s around 20 FPS, which is responsive enough to play with even if it’s not perfectly smooth. The biggest wins:

Persistent thread pool. The original code created and joined OS threads for every parallel_for call — roughly 45 thread creation cycles per frame. Replacing this with a persistent pool (workers sleep on a condition variable, wake on dispatch) roughly halved simulation time. This single change was worth more than all other optimizations combined.

Gradient caching. The spiky kernel gradient \nabla W_{ij} is computed for every neighbor pair during density evaluation. Since the same gradient is needed in both the density correction and divergence correction loops (up to 8 iterations each), caching it once and reusing it across all iterations eliminates millions of redundant sqrt calls per frame.

ZB evaluation skip. The Zhu-Bridson level set is expensive (~17ms). With temporal blending at 0.8, skipping every other frame is nearly invisible but saves half that cost.

Splashes: where particle-only methods hit the wall

The simulation handles sloshing beautifully — waves, vortices, the water climbing up the bunny’s ears and falling back. But splashes are where we hit the fundamental limit of purely Lagrangian SPH.

In a real splash, droplets detach from the bulk, fly through the air, and either rejoin the surface or evaporate. In SPH, a particle that leaves the bulk suddenly has very few neighbors. Its density estimate drops, the pressure solver sees this as “low density” and either ignores the particle or tries to push it back. There’s no clean mechanism for a particle to leave the surface intentionally.

We implemented heuristics: particles with fewer than 12 neighbors are classified as “spray” and exempted from pressure correction, letting them follow ballistic trajectories. Surface particles get an outward velocity boost on impact. This produces visible spray, but it’s a hack — the particles don’t know they’ve detached, they’re just temporarily uncorrected.

We also tried a dual-particle ISPH approach (Liu et al. 2024) that uses virtual particles in the air region for proper bilateral pressure correction. But it finally turned out to require an auxiliary grid — defeating the purpose for our use case. Truly convincing splashes without a grid would need explicit particle splitting and merging, which is a different class of method.

What Claude Code contributed

This project was built entirely with Claude Code as implementation partner. The workflow was iterative: I’d describe what I wanted (a solver, a renderer, an optimization), discuss the approach, and Claude Code would write the implementation. When things broke — and they broke often — we’d debug together, with Claude Code reading error output, analyzing numerical behavior, and proposing fixes.

Claude Code was particularly useful for:

• Implementation of known algorithms: Given a paper reference, Claude Code could produce a working DFSPH solver, Marching Cubes implementation, or caustic renderer without much hand-holding.
• Performance optimization: Essentially all the optimizations — the persistent thread pool, gradient caching, ZB frame skipping — came from Claude Code. Given a profiling picture of where time was being spent, it could identify the bottleneck, propose a fix, and implement it. The 8 → 20 FPS improvement was almost entirely Claude Code’s work.
• Research: When the density field failed near walls, Claude Code helped survey alternatives and identified the Zhu-Bridson approach.

Key takeaways

Auto-compute what you can. The rest density bug would have been nearly impossible to diagnose from visual symptoms alone. Making the solver self-calibrating eliminated an entire class of configuration errors.

Surface reconstruction needs its own attention. A working solver doesn’t give you a nice-looking surface for free. Getting the water to meet container walls cleanly and not jitter frame-to-frame required a complete change of approach (density field to Zhu-Bridson).

Know when you’ve hit the method’s limits. Purely particle-based SPH can simulate sloshing, waves, and gentle splashing convincingly. For dramatic spray and droplet formation, you need either grid-hybrid methods or explicit particle management. Recognizing this boundary saved us from chasing diminishing returns.

Try it yourself

The full source code is available at https://github.com/ai-sci-computing/sloshing-bunny. Build with CMake, run, and shake the bunny.

cmake -B build
cmake --build build
./build/sloshing-bunny

Controls:

• Left drag: Rotate the container (arcball)
• Right drag: Translate (creates inertial sloshing)
• Scroll: Zoom
• Arrow keys: Orbit the camera

References

• Bender, J. & Koschier, D. (2015). Divergence-Free Smoothed Particle Hydrodynamics. Symposium on Computer Animation.
• Macklin, M. & Muller, M. (2013). Position Based Fluids. ACM Trans. Graphics (SIGGRAPH).
• Zhu, Y. & Bridson, R. (2005). Animating Sand as a Fluid. ACM Trans. Graphics (SIGGRAPH).
• Muller, M., Charypar, D. & Gross, M. (2003). Particle-Based Fluid Simulation for Interactive Applications. Symposium on Computer Animation.
• Wyman, C. (2005). An Approximate Image-Space Approach for Interactive Refraction. GPU Gems 2, Chapter 2.