Fluid Simulation

The Fluid Grid

Fluid simulation on the GPU typically uses an Eulerian approach: instead of tracking individual fluid particles, we divide space into a fixed grid and track properties at each cell. Each cell stores velocity (how fast and in what direction fluid is moving), density (how much "stuff" is there), and potentially other quantities like temperature or color.

This grid-based approach maps naturally to GPU texture operations. A 2D fluid simulation stores two textures—one for velocity (two components per pixel) and one for density. Each compute pass reads from one texture and writes to another, iteratively evolving the simulation.

Click and drag to inject fluid and impart velocity. The fluid responds, swirls, and gradually dissipates. This behavior emerges from solving the Navier-Stokes equations on the grid.

The Navier-Stokes Equations

The Navier-Stokes equations describe how fluids move. They look intimidating in mathematical form, but their meaning is physical and intuitive:

This equation says that the change in velocity over time depends on four things. The advection term, $-(\vec{u} \cdot \nabla)\vec{u}$, describes fluid carrying itself along—if there's a current flowing right, that current carries the velocity rightward too. The pressure term, $-\frac{1}{\rho}\nabla p$, pushes fluid from high pressure to low pressure. The viscosity term, $\nu \nabla^2 \vec{u}$, makes neighboring velocities blend together, creating thickness and smoothness. External forces, $\vec{f}$, add gravity, user input, or any other influence.

We also enforce incompressibility: $\nabla \cdot \vec{u} = 0$. This constraint says that fluid cannot compress or expand—what flows into a region must flow out. Without this constraint, fluid would bunch up in some areas and thin out in others, which does not match how water or most liquids behave.

Rather than solving the full equation analytically, we split it into separate steps that we can solve iteratively.

Advection

Advection moves quantities through the velocity field. If the velocity at a point is moving rightward, the density (or velocity itself) at that point should come from whatever was to the left a moment ago.

The standard approach is semi-Lagrangian advection with backward tracing. For each cell, we ask: "Where did the fluid at this cell come from?" We trace backward along the velocity vector, find the source position, and interpolate the value from neighboring cells at that source.

The backtrace finds where fluid came from. We then bilinearly interpolate between the four neighboring cells at the source location. This method is unconditionally stable—it will not explode regardless of timestep—though large timesteps may lose detail.

@compute @workgroup_size(8, 8)
fn advect(@builtin(global_invocation_id) id: vec3<u32>) {
  let pos = vec2<f32>(id.xy) + 0.5;
  let vel = textureLoad(velocity, id.xy, 0).xy;
  
  // Trace backward
  let source = pos - vel * dt;
  
  // Bilinear interpolation at source
  let new_val = textureSampleLevel(field, sampler, source / resolution, 0);
  
  textureStore(output, id.xy, new_val);
}

Self-advection applies the same process to velocity itself: the velocity field is carried along by its own flow. This creates the characteristic swirling behavior of turbulent fluids.

Pressure and Incompressibility

After advection and force application, the velocity field typically has divergence—regions where fluid is compressing or expanding. Real incompressible fluids cannot do this, so we must project the velocity field onto a divergence-free state.

The projection step solves for a pressure field, then subtracts the pressure gradient from velocity. This removes the divergent component, leaving only rotational motion.

The process is:

1. Compute divergence: measure how much fluid is compressing at each cell
2. Solve for pressure using the Poisson equation: $\nabla^2 p = \text{divergence}$
3. Subtract the pressure gradient from velocity: $\vec{u} = \vec{u} - \nabla p$

Step 2 is the expensive part. We use an iterative solver—Jacobi or Gauss-Seidel—that repeatedly averages each cell with its neighbors until the solution converges. More iterations mean better accuracy but cost more compute time.

// Pressure solver iteration
@compute @workgroup_size(8, 8)
fn pressure_iterate(@builtin(global_invocation_id) id: vec3<u32>) {
  let i = vec2<i32>(id.xy);
  
  let div = textureLoad(divergence, i, 0).x;
  
  let pL = textureLoad(pressure, i + vec2(-1, 0), 0).x;
  let pR = textureLoad(pressure, i + vec2(1, 0), 0).x;
  let pB = textureLoad(pressure, i + vec2(0, -1), 0).x;
  let pT = textureLoad(pressure, i + vec2(0, 1), 0).x;
  
  let new_p = (div + pL + pR + pB + pT) / 4.0;
  
  textureStore(pressure_out, i, vec4(new_p, 0.0, 0.0, 0.0));
}

Twenty to forty iterations typically suffice for acceptable visual results. Fewer iterations make the fluid "squishy"; more make it crisper.

A Simple Fluid Solver

The complete solver chains these operations together. Each frame runs multiple compute passes, ping-ponging between buffer pairs to avoid read/write conflicts.

The order matters. We apply forces first, then diffuse (for viscosity), then project to remove divergence. Advection comes after projection so we are moving a divergence-free field. A second projection after advection catches any divergence introduced by numerical errors.

For density, the process is simpler: diffuse (if desired), then advect through the velocity field. Density does not affect velocity in our simple model—it is just carried along for visualization.

Viscosity (diffusion of velocity) uses the same iterative solver as pressure. Higher viscosity makes the fluid thicker, like honey. Zero viscosity gives inviscid flow, which swirls freely without damping.

Rendering the Fluid

The simulation produces fields of numbers. Turning them into visuals is a design choice. Density can render as brightness or color intensity—brighter where there is more fluid. Velocity can render as color mapped by direction, with hue indicating flow angle and saturation indicating speed.

For dye-like effects, track color per cell (RGB instead of scalar density) and advect each channel. For smoke, render density with soft gradients and perhaps alpha blending. For water surfaces, couple the 2D simulation to height-field waves or render a 3D volume.

The rendering pass is typically a fragment shader sampling the simulation textures:

@fragment
fn fs(@location(0) uv: vec2<f32>) -> @location(0) vec4<f32> {
  let density = textureSample(density_tex, sampler, uv).x;
  let velocity = textureSample(velocity_tex, sampler, uv).xy;
  
  // Density as blue intensity
  let color = vec3(0.2, 0.4, 1.0) * density * 2.0;
  
  return vec4(color, 1.0);
}

Boundaries and Extensions

The edges of the grid require special handling. Solid boundaries (walls) reflect velocity: the component perpendicular to the wall is negated. Open boundaries might allow fluid to flow out or wrap around.

The basic solver can be extended in many directions. Vorticity confinement adds small-scale turbulence that would otherwise be lost to numerical diffusion. Temperature can drive buoyancy, making hot fluid rise. Obstacles within the grid require masking—setting velocity to zero inside solids and handling boundary conditions at interfaces.

3D fluid simulation follows the same principles but with three-dimensional grids and correspondingly higher memory and compute costs. The principles—advect, diffuse, project—remain identical; only the neighbor counts and storage requirements change.

Key Takeaways

• Eulerian fluids simulate on a fixed grid, storing velocity and density per cell
• Navier-Stokes equations decompose into advection, pressure, viscosity, and forces
• Semi-Lagrangian advection traces backward through the velocity field
• Pressure projection enforces incompressibility via an iterative Poisson solve
• The solver runs multiple compute passes per frame, ping-ponging between buffers
• Rendering maps simulation fields to colors; density as intensity, velocity as hue