Fluid Dynamics Simulation

physicswasmrustwebglfluid-dynamics

Particle-Based Fluid Simulation

This interactive simulation demonstrates real-time fluid dynamics using the double density relaxation algorithm developed by Clavet, Beaudoin, and Poulin (2005). The entire physics engine is written in Rust and compiled to WebAssembly for maximum performance, achieving 60fps with 3000+ particles.

The simulation features three density classes that stratify naturally: heavy particles (cyan) sink to the bottom, medium particles (blue) occupy the middle, and light particles (purple) float to the top. This demonstrates realistic buoyancy and density-driven stratification as different particle classes separate into distinct horizontal layers.

Loading simulation...

FPS: 0
Frame Time: 0.00ms
Particles: 3000

Click and drag to apply force to the fluid

Algorithm Implementation

The simulation follows the algorithm of Clavet et al. (2005) with extensions for multiple density classes, implementing a three-pass approach each frame:

Pass 1: Viscosity, Forces & Position Update

  • Apply viscosity impulses between approaching particles (Algorithm 5 from paper)
  • Store previous positions for Verlet integration
  • Apply external forces (gravity scaled by density class: 1.2× for heavy, 1.0× for medium, 0.75× for light)
  • Update positions based on velocities
  • Build spatial hashmap for neighbor queries in Pass 2

Pass 2: Double Density Relaxation

  • Freeze positions at start of pass to prevent feedback loops
  • Calculate density (Σg²) and near-density (Σg³) for each particle using smoothing kernels
  • Compute pressure based on deviation from class-specific rest density
  • Apply symmetric pair pressure (P_i + P_j) with mass-weighted displacement splitting
  • Lighter particles move more than heavier ones from same pressure (enables stratification)
  • This creates incompressibility through particle repulsion while preserving momentum

Pass 3: Boundary Containment & Velocity Finalization

  • Constrain particles to circular container with hysteresis threshold (0.5px dead zone)
  • Snap particles slightly inside boundary to prevent sticking
  • Reflect outward velocities with 80% energy retention
  • Calculate final velocities from position delta (Verlet integration)

Critical Implementation Details

Immediate Neighbor Updates: The algorithm processes each particle sequentially, modifying neighbor positions immediately rather than accumulating all changes. This creates proper bidirectional forces where each pair of particles pushes each other based on their individual pressures.

Mass-Weighted Forces: Displacement magnitudes are split proportionally by inverse mass (1/rest_density), allowing lighter particles to move more from the same pressure differential. This enables realistic stratification where heavy particles accumulate at the bottom.

Frozen Reference Positions: Density calculations in Pass 2 use frozen position copies to prevent cascading modifications during relaxation, ensuring algorithmic correctness and numerical stability.

Viscosity Timing: Viscosity impulses are applied in Pass 1 using the previous frame’s spatial hashmap, before clearing for new positions. This follows Algorithm 5 from the paper and prevents mixing that would counteract stratification.

Spatial Hashing: An efficient O(n) neighbor-finding system using spatial grid partitioning with domain offset to handle negative coordinates correctly, avoiding O(n²) all-pairs checks.

Performance Optimizations

Rust/WASM Pipeline:

  • Zero-copy memory access via shared TypedArrays
  • Structure of Arrays (SoA) layout for cache coherency
  • SIMD-friendly contiguous data (all X positions together, all Y positions together)
  • Compiled to near-native performance with LTO and optimization level 3

WebGL 2 Rendering:

  • GPU-instanced particle rendering (single draw call for all particles)
  • Alpha blending for consistent opacity (prevents brightness stacking)
  • Soft-edged particles with smooth falloff for blob-like appearance
  • 60fps with 3000+ particles, stable up to 4000
  • Fermat spiral (sunflower pattern) initialization for uniform distribution

Efficient Algorithms:

  • Fast approximate normalization (avoiding sqrt when possible)
  • Spatial grid with O(1) cell lookups
  • Early-exit checks for zero forces
  • Comprehensive bounds checking to prevent crashes

Physics Parameters

The simulation uses constants derived from the Clavet paper and extensive tuning:

  • Interaction Radius: 30.0 pixels
  • Time Step: 0.0166 seconds (60fps)
  • Stiffness: 0.004 / dt² ≈ 14.5 (scaled for small timesteps per Clavet Algorithm 2)
  • Near Stiffness: 0.010 / dt² ≈ 36.3
  • Rest Densities: Heavy=12.0, Medium=10.0, Light=8.0 (class-specific for stratification)
  • Gravity: 200.0 pixels/s² scaled by class (1.2×, 1.0×, 0.75×)
  • Container Radius: 400.0 pixels
  • Boundary Damping: 0.8 (retains 80% of velocity on collision)
  • Boundary Threshold: 0.5 pixels (dead zone prevents jitter)
  • Viscosity: σ=0.02, β=0.05 (reduced to allow density separation)

Interaction Controls

Click and drag anywhere on the canvas to apply attracting or repelling forces to the fluid. Particles respond dynamically with realistic fluid behaviors:

  • Incompressibility: Particles maintain target density through pressure forces
  • Surface Tension: Near-density forces create cohesive blobs
  • Momentum Conservation: Forces obey Newton’s third law
  • Density Stratification: Different particle classes naturally separate
  • Viscosity: Velocity damping creates realistic fluid friction

Technical Stack

  • Rust - Core physics simulation with comprehensive bounds checking
  • WebAssembly - Compiled from Rust for near-native browser performance
  • WebGL 2 - GPU-accelerated instanced rendering
  • React - UI controls and state management
  • TypeScript - Type-safe integration layer

Development Insights

This implementation went through extensive debugging to achieve physically accurate stratification. Key challenges solved:

  1. Stiffness Scaling: Clavet’s paper uses dt=1, but at dt≈0.0166 (60fps), stiffness must scale by 1/dt² (~3600×) to maintain equivalent pressure forces. Missing this scaling caused particles to act like gas instead of liquid.

  2. Mass-Weighted Displacement: Equal 50/50 force splitting prevents stratification because heavy and light particles move equally. Proper mass-weighted splitting (displacement ∝ inverse mass) is critical for density classes to separate.

  3. Viscosity Hashmap Bug: Applying viscosity with an empty hashmap (cleared too early) resulted in zero viscosity forces. Correct ordering: apply viscosity → clear hashmap → update positions → rebuild hashmap.

  4. Boundary Jitter: Particles oscillating across boundary threshold created visible jitter. Solution: hysteresis dead zone (0.5px threshold), snap close to boundary (399.8px), and high velocity retention (80%) for smooth behavior.

  5. Rendering Mode: Additive blending (gl.ONE) caused brightness stacking where particles overlapped, masking the true density distribution. Alpha blending (gl.ONE_MINUS_SRC_ALPHA) maintains consistent color.

  6. Particle Initialization: Uninitialized old_x/old_y positions caused wild velocities on first frame (Verlet integration bug). All initialization methods must set old positions equal to current positions.

  7. Coordinate System Confusion: Renderer Y-axis flip meant gravity direction was inverted. Positive Y in simulation → negative Y on screen after flip.

References

Primary Source:

Clavet, S., Beaudoin, P., & Poulin, P. (2005). Particle-based Viscoelastic Fluid Simulation.
Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation.

Implementation Reference:

The algorithm follows “Algorithm 2: Double density relaxation” and “Algorithm 5: Viscosity impulses” from the original paper, with careful attention to the specific update ordering and force calculations described.

Architecture

This simulation uses a unified monorepo architecture where all source code (Rust physics engine, TypeScript renderer, and React UI wrapper) lives in the same repository. This eliminates synchronization issues between separate codebases and ensures the physics implementation and visual rendering stay perfectly aligned.

This Rust/WASM implementation achieves significant performance improvements over JavaScript implementations while maintaining exact algorithmic fidelity to the reference paper, with additional enhancements for multi-density stratification.