What is Fluid Simulation?

Fluid simulation is the computational modeling of how liquids and gases move. From the swirling patterns in your coffee to the formation of hurricanes, fluids follow mathematical laws that we can simulate using computers.

The challenge? Fluids are continuous - they have properties at every point in space. Computers are discrete - they can only work with finite numbers. Our task is to approximate continuous fluid behavior using discrete computations.

The Navier-Stokes Equations

At the heart of fluid simulation lie the Navier-Stokes equations, formulated in the 19th century. These partial differential equations describe how fluid velocity changes over time.

Operator Splitting: Divide and Conquer

Solving the full Navier-Stokes equations directly is computationally expensive. Instead, we use operator splitting to break the problem into simpler steps that we solve sequentially.

Each step solves a simpler equation that, when combined, approximates the full Navier-Stokes behavior.

Advection: Following the Flow

Advection describes how quantities (velocity, color, temperature) are transported by the fluid flow. The advection equation is:

$$\frac{\partial q}{\partial t} + \mathbf{u} \cdot \nabla q = 0$$

Where \(q\) is any quantity we're transporting (like dye color or velocity itself).

Semi-Lagrangian Advection

We use the Semi-Lagrangian method, which is unconditionally stable. The idea:

1. For each grid point, trace backwards along the velocity field
2. Find where this particle came from at the previous timestep
3. Sample the quantity at that location (with interpolation)

// Backtrace particle position
vec2 coord = vUv - dt * texture2D(uVelocity, vUv).xy * texelSize;

// Sample quantity at previous position
vec4 result = texture2D(uSource, coord);

// Apply dissipation
float decay = 1.0 + dissipation * dt;
gl_FragColor = result / decay;

Pressure Projection: Enforcing Physics

After advection and forces, our velocity field likely has divergence (\(\nabla \cdot \mathbf{u} \neq 0\)), violating incompressibility. Pressure projection fixes this.

The Method

We use the Helmholtz-Hodge decomposition, which says any vector field can be split into divergence-free and curl-free components. We want the divergence-free part.

Jacobi Iteration

The Poisson equation is solved iteratively using the Jacobi method:

$$p_{i,j}^{k+1} = \frac{p_{i+1,j}^k + p_{i-1,j}^k + p_{i,j+1}^k + p_{i,j-1}^k - \text{div}_{i,j}}{4}$$

float L = texture2D(uPressure, vL).x;
float R = texture2D(uPressure, vR).x;
float T = texture2D(uPressure, vT).x;
float B = texture2D(uPressure, vB).x;
float divergence = texture2D(uDivergence, vUv).x;

// One Jacobi iteration
float pressure = (L + R + B + T - divergence) * 0.25;
gl_FragColor = vec4(pressure, 0.0, 0.0, 1.0);

We perform 20-40 iterations per frame, progressively reducing divergence. More iterations = more accurate.

GPU Implementation: Massive Parallelism

Modern GPUs can perform thousands of computations simultaneously. We exploit this by representing our simulation grid as textures and computations as fragment shaders.

Why GPU?

• Parallelism: 16,384 pixels = 16,384 simultaneous calculations
• Memory Locality: Textures naturally represent 2D grids
• Hardware Optimization: GPUs are built for this workload

Ping-Pong Rendering

Shaders can't read and write to the same texture. We use ping-pong buffering: two textures, alternating between read and write.

// Read from one FBO, write to another
shader.bind();
gl.uniform1i(uniforms.uSource, velocity.read.attach(0));
fboManager.blit(velocity.write);

// Swap for next frame
velocity.swap(); // read ↔ write

Resolution Strategy

We run physics at low resolution (128×128) and visuals at high resolution (1024×1024). This decoupling gives us both performance and quality.

Visual Effects: Making It Beautiful

Raw simulation output looks flat. We add post-processing to enhance aesthetics without affecting physics.

Bloom Effect

Bloom creates a glow around bright areas, mimicking how cameras and human eyes see intense light.

1. Prefilter: Extract pixels above brightness threshold
2. Blur: Apply Gaussian blur at multiple resolutions
3. Composite: Add blurred glow to original image

Sunrays (God Rays)

Radial blur creates light streaks from bright sources, like sunlight through clouds.

// Sample along ray from center
for (int i = 0; i < 16; i++) {
    coord -= dir * density;
    color += texture2D(uTexture, coord).a * illuminationDecay;
    illuminationDecay *= decay;
}

Shading

We compute surface normals from the dye density gradient, then apply simple lighting for depth perception.

Obstacles: Solid Boundary Conditions

Real fluids don't just flow freely - they interact with solid surfaces. Our simulation includes an obstacle (the "M" shape you see) that demonstrates how fluids behave around boundaries.

The Challenge

Solids impose two key physical constraints:

• No-Penetration: Fluid cannot pass through walls
• No-Slip: Viscous fluid at a wall has zero velocity

Mathematically, at a solid boundary:

$$\mathbf{u}|_{\text{wall}} = 0$$

Representing Obstacles

We use a binary obstacle mask texture - each pixel stores whether that location is solid (1) or fluid (0). This mask is sampled in every shader to modify behavior near boundaries.

// Create obstacle texture from rectangles
const mask = new Uint8Array(width * height);

for (const obstacle of obstacles) {
    for (let y = 0; y < height; y++) {
        for (let x = 0; x < width; x++) {
            if (isInsideRectangle(x, y, obstacle)) {
                mask[y * width + x] = 255;  // Mark as obstacle
            }
        }
    }
}

Boundary Conditions in Each Solver Step

1. Advection

When back-tracing particle paths, we check if the traced position lands inside an obstacle. If so, we clamp it to the boundary - particles can't advect through walls.

vec2 tracedPos = vUv - dt * velocity * texelSize;

// Check if traced position hits obstacle
if (texture2D(uObstacles, tracedPos).r > 0.5) {
    tracedPos = vUv;  // Stay at current position
}

2. Divergence Computation

The most subtle part. When computing \(\nabla \cdot \mathbf{u}\), we sample neighboring velocities. If a neighbor is an obstacle, we apply the no-slip reflection:

$$u_{\text{obstacle}} = -u_{\text{current}}$$

This reflection ensures the average velocity at the boundary is zero.

float sampleVelocity(vec2 coords, float currentComponent) {
    // Domain boundaries
    if (coords.x < 0.0 || coords.x > 1.0 || 
        coords.y < 0.0 || coords.y > 1.0) {
        return -currentComponent;  // Reflect at domain edge
    }
    
    // Check for obstacle
    if (texture2D(uObstacles, coords).r > 0.5) {
        return -currentComponent;  // No-slip: reflect
    }
    
    // Normal fluid cell
    return texture2D(uVelocity, coords);
}

3. Pressure Solve

Obstacles use Neumann boundary conditions: \(\frac{\partial p}{\partial n} = 0\). The pressure gradient normal to the wall is zero.

In practice, when sampling pressure neighbors during Jacobi iteration, if a neighbor is an obstacle, we use the current cell's pressure (no gradient).

float samplePressure(vec2 coords) {
    if (texture2D(uObstacles, coords).r > 0.5) {
        // Neumann BC: return current cell pressure
        return texture2D(uPressure, vUv).r;
    }
    return texture2D(uPressure, coords).r;
}

4. Gradient Subtraction

After solving for pressure, we subtract \(\nabla p\) from velocity. In obstacle cells, we skip this and simply enforce zero velocity.

vec2 u = texture2D(uVelocity, vUv).xy;

if (texture2D(uObstacles, vUv).r > 0.5) {
    gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);  // No velocity in obstacles
} else {
    vec2 gradient = vec2(pR - pL, pT - pB) * 0.5;
    gl_FragColor = vec4(u - gradient, 0.0, 1.0);
}

Visualizing Obstacles

The obstacle itself is rendered in the final display pass by sampling the mask texture and blending a distinct color (bright cyan) where obstacles exist.

vec3 color = texture2D(uDye, vUv).rgb;

#ifdef SHOW_OBSTACLES
    float obstacle = texture2D(uObstacles, vUv).r;
    if (obstacle > 0.5) {
        color = uObstacleColor;  // Override with obstacle color
    }
#endif

gl_FragColor = vec4(color, 1.0);

Challenges & Solutions

Numerical Stability

Problem: Large timesteps cause explosions.

Solution: Semi-Lagrangian advection is unconditionally stable. We clamp timesteps to 16ms max.

Numerical Dissipation

Problem: Interpolation and iteration smooth out details.

Solution: Vorticity confinement adds energy back to small vortices.

Browser Compatibility

Problem: Different browsers support different WebGL features.

Solution: Detect capabilities, fallback from WebGL2 to WebGL1, handle missing extensions gracefully.

Performance Optimization

• Half-Float Textures: 16-bit precision instead of 32-bit
• Resolution Scaling: Adapt to device capabilities
• Iteration Budgets: Fewer pressure iterations on mobile
• Effect Toggles: Disable bloom/sunrays on low-end devices

Bug Fixes & Improvements (April 2026)

After a thorough code review, several bugs and improvements were identified across the simulation pipeline. Here's a summary of the fixes applied:

Gradient Subtraction Scaling Fix

The pressure projection step had an inconsistency: the divergence shader correctly used central differences with a 0.5 factor (div = 0.5 * (R - L + T - B)), but the gradient subtraction shader omitted it. This made the pressure correction twice as strong as intended, causing over-correction artifacts. The fix:

// Before (incorrect)
velocity.xy -= vec2(R - L, T - B);

// After (correct - matches divergence scaling)
velocity.xy -= 0.5 * vec2(R - L, T - B);

Dissipation Formula Upgrade

The advection shader used a first-order Taylor approximation for dissipation: result / (1.0 + dissipation * dt). This is inaccurate for larger dissipation values and doesn't model true exponential decay. Replaced with the physically correct formula:

// Before (approximation)
float decay = 1.0 + dissipation * dt;
gl_FragColor = result / decay;

// After (exact exponential decay)
gl_FragColor = result * exp(-dissipation * dt);

Pressure Solver Convergence

The Jacobi pressure solver defaulted to 20 iterations, which is insufficient for a 256x256 grid to reach adequate convergence. Residual divergence was leaving compressibility artifacts in the velocity field. Increased the default to 40 iterations, which roughly halves the residual error without excessive GPU cost.

Timestep Double-Clamping

Delta time was being clamped to 16ms in both the render loop (main.js) and inside SimulationManager.update(). The redundant clamp was removed. At low frame rates (e.g., 30 FPS), the double-clamp was causing the simulation to appear to run at half real-time speed.

Touch Pointer Memory Leak

The PointerManager tracked multi-touch pointers by adding entries to an array on touchstart, but never removed them on touchend — only setting down = false. Over extended touch-heavy sessions, the array grew unboundedly. Stale pointers are now cleaned up when touches end.

Other Fixes

• Duplicate method removed: _loadShaderFile was defined twice in SimulationManager due to a copy-paste error.
• FBO buffer references: The full-screen quad's vertex and index buffers are now stored as quadVBO/quadIBO so they can be properly cleaned up on WebGL context loss.
• Color loop cleanup: The pointer color transition's requestAnimationFrame loop now stores its ID and exposes a destroy() method, preventing it from running indefinitely after the simulation is torn down.

Conclusion

Real-time fluid simulation combines classical physics, numerical methods, and modern GPU programming. By breaking the Navier-Stokes equations into manageable pieces and leveraging parallel hardware, we achieve beautiful, interactive simulations in the browser.

The techniques explored here extend beyond fluid simulation - they're fundamental to real-time physics, weather modeling, game engines, and scientific visualization.

Further Reading

• Stam, Jos. "Stable Fluids." SIGGRAPH 1999.
• Harris, Mark J. "Fast Fluid Dynamics Simulation on the GPU." GPU Gems Chapter 38, NVIDIA, 2004.
• Bridson, Robert. "Fluid Simulation for Computer Graphics." A K Peters/CRC Press, 2008.
• Fedkiw, Ronald, et al. "Visual Simulation of Smoke." SIGGRAPH 2001.

Acknowledgments

This implementation is heavily inspired by Pavel Dobryakov's excellent WebGL Fluid Simulation . His work provided invaluable insights into GPU-based fluid dynamics and modern WebGL techniques.

Big thanks to Pavel for creating such an elegant and well-engineered implementation that serves as an outstanding reference for the fluid simulation community.

Further Reading

• Stam, Jos. "Stable Fluids." SIGGRAPH 1999.
• Harris, Mark J. "Fast Fluid Dynamics Simulation on the GPU." GPU Gems Chapter 38, NVIDIA, 2004.
• Bridson, Robert. "Fluid Simulation for Computer Graphics." A K Peters/CRC Press, 2008.
• Fedkiw, Ronald, et al. "Visual Simulation of Smoke." SIGGRAPH 2001.