realtime-cloth-fluid

How it is built

Architecture overview

A real-time cloth and fluid simulator advances thousands to millions of points by a small fixed time step, several times per rendered frame, and must never blow up. The one idea that makes this possible is projection instead of force integration: cloth and soft bodies are particles whose predicted positions are pushed back onto constraints (distance, bending, volume), and a liquid is made incompressible by projecting its velocity onto a divergence-free field or its density onto a target. A large step then gives an inaccurate but bounded result, not an explosion. Around that core sit the machinery of any particle system: a neighbour search so particles find each other without an all-pairs loop, collision detection against the world, a linear solver for pressure, and a surface extractor that turns particles or a field into a mesh.

Parts:

• Time stepping — fixed-step loop and the integrator that predicts positions — Semi-implicit Euler, built on fixed-timestep-loop (no page yet)
• Cloth constraint solver — projects predicted positions onto constraints, in parallel batches — position-based-dynamics (no page yet), graph-coloring (no page yet), built on BFS and DFS and d-ary heap
• Neighbour search — finds the particles near each particle every step — Spatial hashing, built on Fast hash functions; Radix sort; parallel-prefix-sum (no page yet)
• Collision — broad phase against meshes, and time of impact so thin cloth is not tunnelled through — Bounding volume hierarchy; continuous-collision-detection (no page yet), built on GJK and EPA, Newton root finding and Sweep and prune
• Particle fluid — liquid as free particles with density-based pressure — smoothed-particle-hydrodynamics (no page yet), built on neighbor-lists (no page yet) and Verlet and leapfrog
• Grid fluid — liquid or smoke as a velocity field on a grid, particles carrying the liquid — flip-pic-fluid (no page yet), built on PCG and xoshiro generators; semi-lagrangian-advection (no page yet), built on Runge-Kutta 4 and spline-interpolation (no page yet); pressure-projection-method (no page yet), built on Finite differences and Incomplete factorization preconditioners
• Pressure solver — the linear system of the grid fluid — Conjugate gradient, Multigrid
• Surface — turns particles and fields into a renderable surface — level-set-methods (no page yet), built on flux-limiters-upwinding (no page yet); marching-cubes (no page yet)

 user input / animation ──► fixed-step loop (n substeps per frame)
                                  │
              ┌───────────────────┼───────────────────────────┐
              ▼                   ▼                           ▼
        cloth / soft body    particle fluid (SPH/PBF)    grid fluid (FLIP)
        predict ─► project   hash grid ─► density ─►     particles ─► grid ─► advect
        constraints (colors) pressure ─► update           ─► pressure (CG/MG) ─► particles
              │                   │                           │
              └──── collision: BVH broad phase, swept test ◄──┘
                                  │
                      level set / marching cubes ─► mesh ─► render

The three solvers share the time loop, the neighbour search and the collision code. A game typically ships one or two of them, not all three.

How the parts fit together

Semi-implicit Euler

Semi-implicit Euler updates velocity from forces first, then position from the new velocity: explicit-Euler cost, but symplectic, so it is the default game-physics integrator. Here it is the prediction step that PBD, PBF and FLIP all follow with a correction. Contract with the loop: fixed-timestep-loop (no page yet) hands every solver the same constant dt from an accumulator and may split a frame into substeps. Verlet and leapfrog is time-reversible and second order and serves the particle fluid; the cloth path keeps semi-implicit Euler because constraint projection overwrites accuracy anyway.

Position based dynamics

position-based-dynamics (no page yet) is the cloth and soft-body solver. After the prediction each constraint moves the predicted positions of its particles to satisfy it, then velocities are recomputed from the position change. Plain PBD ties stiffness to iteration count and step; XPBD adds a compliance and a per-constraint multiplier so stiffness becomes a material parameter independent of both (Macklin, Müller, Chentanez, 2016). The alternative is the force-and-impulse route of Sequential impulse solver, better for rigid stacks and joints; for a flag or soft body the position-level solver is simpler and more forgiving. Contract: it takes predicted positions and masses from the integrator, contact constraints from collision, and returns corrected positions.

Graph colouring

graph-coloring (no page yet) makes the sequential Gauss-Seidel sweep parallel. Constraints that share no particle can be solved at the same time without a race, so colour the constraint graph (a greedy or DSatur colouring; BFS and DFS supply the traversal and DSatur needs a priority queue such as a d-ary heap) and launch one parallel kernel per colour. XPBD's paper formulates both Gauss-Seidel and Jacobi updates and shows both converge linearly; colouring keeps the Gauss-Seidel form. Vivace (Fratarcangeli, Tibaldo, Pellacini, 2016) is a related route, a parallel randomized Gauss-Seidel method. Contract: the solver receives a list of batches; inside a batch no order is assumed. The colouring is rebuilt only when constraints appear.

Spatial hashing

Spatial hashing keeps a uniform grid in a hash table: a particle goes into the cell of its position and a query reads the surrounding cells. With cell size equal to the interaction radius the search is O(n) for similar-sized particles, which is the case for SPH particles and for cloth vertices in self-collision. A tree adapts to varying sizes but costs more per query and per rebuild. Fast hash functions map the integer cell coordinates to a slot; a table smaller than the grid only adds false candidates. Contract: rebuilt each step from current positions; it returns candidates, and consumers filter by true distance.

Radix sort

Radix sort builds the same grid on a GPU. NVIDIA's CUDA particle whitepaper computes a cell hash per particle, sorts by it with a radix sort, and then finds the start of each cell in the sorted array (Green, 2010). Compared with its atomic-counter grid, sorting needs no fixed maximum per cell and no contended atomics, and it places spatially close particles close in memory. Contract: the output is a permutation plus a cell-start array; state arrays are reordered, or read through the permutation, before the density or collision pass.

Parallel prefix sum

parallel-prefix-sum (no page yet) turns per-cell counts into start offsets, and compacts lists (for example only active contacts) without a sequential loop. The same whitepaper describes its counting-sort grid build as a prefix sum that resembles one pass of a radix sort. It is shared infrastructure of the neighbour search, the cloth contact list and the particle emitter. Contract: plain arrays in, plain arrays out.

Bounding volume hierarchy

Bounding volume hierarchy is the broad phase against meshes (a character, a table, a level). It prunes the triangles that cannot touch a cloth vertex to a logarithmic query and copes with objects of very different sizes, which a uniform grid does not; after each step the boxes are refitted instead of rebuilt. Contract: collision asks for pairs (cloth element, world triangle) whose boxes overlap, and hands them to the narrow phase.

Continuous collision detection

continuous-collision-detection (no page yet) finds when during the step two moving elements first touch. Cloth is thin, so a fast vertex can pass through a triangle between steps and an end-of-step test never notices; cloth collision work therefore uses swept tests (Bridson, Fedkiw, Anderson, 2002). The time of impact is a root-finding problem (Newton root finding); convex shapes can use GJK and EPA; candidate pairs are first reduced by Sweep and prune or the hierarchy above. Many real-time systems cap the motion per substep and add a margin, and use continuous tests only for fast objects. Contract: it reports time and normal, which the cloth solver turns into a position constraint.

Smoothed particle hydrodynamics

smoothed-particle-hydrodynamics (no page yet) models a liquid as free particles; density, pressure and viscosity are sums over neighbours weighted by a kernel. There is no grid, so splashes, free surfaces and coupling with cloth particles come cheaply. Simple SPH needs small steps to stay stiff, so libraries use incompressible variants: SPlisHSPlasH lists WCSPH, PCISPH, PBF, IISPH, DFSPH and PF; Position Based Fluids (Macklin, Müller, 2013) puts the density constraint into the cloth solver's framework. The alternative, FLIP, behaves better for large water volumes but needs a grid and a pressure solve. Contract: candidates come from the hash grid; neighbor-lists (no page yet) cache them across substeps; Verlet and leapfrog moves the particles.

FLIP and PIC

flip-pic-fluid (no page yet) is the hybrid: particles carry mass and velocity; each step they transfer velocity to a grid, the grid solves pressure, and the velocity change goes back to the particles (FLIP) or the whole new velocity (PIC). Blending the two trades the diffusion of PIC against the noise of FLIP; APIC also transfers a local affine velocity. Zhu and Bridson (2005) introduced it to graphics for sand and water. Choose it when volume and depth should look right, SPH for small splashy liquids. Jittered seeding per cell uses a PCG or xoshiro generator: fast, and a replay gives the same liquid. Contract: the two transfers are the interface to the grid stages; cells are marked liquid or air from the particles.

Semi-Lagrangian advection

semi-lagrangian-advection (no page yet) moves velocity, smoke density or temperature through the grid by tracing each sample backwards along the velocity and interpolating there. Stam's Stable Fluids (1999) built real-time fluids on it because any dt is stable. The backtrace uses Runge-Kutta 4 or a midpoint step, the lookup spline-interpolation (no page yet), linear for speed, cubic for less blur. The price is numerical diffusion; smoke adds vorticity confinement (Fedkiw, Stam, Jensen, 2001), and in FLIP the particles carry the detail. Contract: velocity field in, advected field out, which must then be made incompressible.

Pressure projection

pressure-projection-method (no page yet) is Chorin's projection (1968): solve a Poisson equation for the pressure whose gradient removes the divergence of the velocity, then subtract that gradient. Finite differences on a staggered (MAC) grid keep divergence and gradient consistent, and solid and air cells enter as boundary conditions. Taichi's stable_fluid.py example does exactly this after semi-Lagrangian advection, with p_jacobi_iters = 500 by default (40 for a quicker, less accurate result in its comment) or a sparse-matrix solve. Contract: it hands a sparse symmetric system to the pressure solver, which may use a preconditioner from Incomplete factorization preconditioners.

Conjugate gradient

Conjugate gradient solves that system matrix-free: one stencil application and two dot products per iteration, and about √κ iterations, so the cost climbs with grid resolution. Preconditioning cuts iterations; an incomplete factorization costs a triangular solve that is harder on GPUs. Contract: it needs only a function applying the matrix, so boundary handling stays in the operator. It is the default for CPU games and small grids; Multigrid is its alternative.

Multigrid

Multigrid smooths the high-frequency error on the fine grid and solves the smooth remainder on coarser grids, in work proportional to the number of cells. Large 3D grids, where CG needs many iterations, use it directly or as a preconditioner. The cost is complexity near moving boundaries and free surfaces. Contract: same as CG (system in, pressure out), so the two can sit behind one interface and be compared on the same scene.

Level set methods

level-set-methods (no page yet) represent the liquid surface as the zero set of a signed distance field. The field is advected with the fluid using upwind schemes and slope limiters (flux-limiters-upwinding, no page yet) and a Runge-Kutta step, then reinitialised. It gives normals, curvature and inside tests for surface tension and rendering; for particle liquids a field is built from the particles each frame. Contract: particles or the velocity grid in, a scalar field out; it is a view of the dynamics, not part of them.

Marching cubes

marching-cubes (no page yet) extracts triangles from that field (Lorensen and Cline, 1987): the sign pattern of a cell's eight corners selects a precomputed triangulation. The alternative that skips meshing, rendering particles as screen-space spheres with smoothed depth, is cheap and common in games; choose marching cubes when the surface must be a real mesh (export, collision, shading). Contract: field in, indexed triangles out; normals come from the field gradient (Finite differences).

Data and control flow

One cloth frame (PBD/XPBD):

1. Split the frame into substeps of fixed dt.
2. Semi-implicit Euler adds gravity and forces to velocity and predicts positions.
3. Query the hierarchy (and a hash grid for self-collision) on the predicted positions; contacts become constraints.
4. Per colour, a parallel kernel projects its constraints; a few iterations per substep.
5. Recompute velocities from the position change, apply friction and damping, commit positions.
6. Update the render mesh. Time goes to constraint projection and collision queries.

One FLIP step:

1. Sort particles by cell (radix sort, prefix sum).
2. Transfer velocity to the grid; mark liquid cells.
3. Add gravity; advect other grid quantities.
4. Compute divergence, solve pressure by CG or multigrid, subtract the gradient.
5. Transfer the velocity change back (blended with PIC), advect particles, resolve solid collisions.
6. Per frame, build a level set and extract a mesh, or render particles directly. Time goes to the pressure solve and the transfers.

Scaling

1. Pressure solve (grid fluids) limits first. Levers: a preconditioner, a warm start from the last frame, multigrid, a coarser grid with particles for detail.
2. Neighbour search dominates SPH and self-collision. Levers: cell size equal to the radius, sorted particle layout, neighbour lists reused across substeps.
3. Constraint count bounds cloth. Levers: more substeps with fewer iterations, colouring for parallel kernels, a coarse simulation mesh driving a detailed render mesh.
4. Collision grows with scene complexity. Levers: refitted hierarchy, margins, swept tests only for fast objects.
5. GPU memory bandwidth: structure-of-arrays state, sorted particle order, 32-bit floats.

Failure modes

Real implementations

The alternative: offline VFX solvers use fine grids and higher-order schemes and are not bound to a frame budget. This design trades accuracy for bounded cost per frame; when results must be physically accurate use a validated finite-element or CFD solver.

Sources

• M. Müller, B. Heidelberger, M. Hennix, J. Ratcliff, Position based dynamics, 2007. DOI
• M. Macklin, M. Müller, N. Chentanez, XPBD, MIG 2016. paper
• M. Macklin, M. Müller, Position based fluids, 2013. DOI
• J. Stam, Stable fluids, 1999. DOI
• Y. Zhu, R. Bridson, Animating sand as a fluid, 2005. DOI
• A. J. Chorin, Numerical solution of the Navier-Stokes equations, 1968. DOI
• R. Fedkiw, J. Stam, H. W. Jensen, Visual simulation of smoke, 2001. DOI
• R. Bridson, R. Fedkiw, J. Anderson, Robust treatment of collisions, contact and friction for cloth animation, 2002. DOI
• M. Fratarcangeli, V. Tibaldo, F. Pellacini, Vivace, 2016. DOI
• W. Lorensen, H. Cline, Marching cubes, 1987. DOI
• S. Green, Particle Simulation using CUDA, NVIDIA, 2010. PDF
• Taichi example stable_fluid.py