Learn & Review: Coding Adventure: Simulating Fluids

Coding Adventures: Fluid Simulations with Smoothed Particle Hydrodynamics

This summary outlines the process of building a fluid simulation using Smoothed Particle Hydrodynamics (SPH), detailing the challenges encountered and the solutions developed.

1. Initial Setup and Basic Physics

• Core Concept: Representing a fluid as a collection of particles that move in response to forces.
• Initial Particle: A single circle representing a tiny bit of fluid.
• Forces Introduced:
  • Gravity: Accelerates particles downwards.
  • Collisions: Particles bounce off boundaries of a defined box.
  • Collision Damping: Controls energy loss during bounces, causing particles to eventually stop.

2. Scaling to Multiple Particles and the Overlapping Problem

• Multiple Particles: The simulation is scaled up by using arrays to store positions and velocities for many particles.
• Particle Arrangement: Particles are initially arranged in a grid for better distribution.
• The Problem: Particles collapse on top of each other due to a lack of repulsive forces.

3. Introducing Smoothed Particle Hydrodynamics (SPH)

• Goal: To simulate fluid behavior, specifically addressing the overlapping particle issue.
• Core Idea of SPH: To approximate a continuous field (like density) from discrete particles by "blurring" or "smoothing" their influence over a region. This technique was originally used in astrophysics.
• Smoothing Function:
  • Defines a "circle of influence" around each particle.
  • Influence is maximum at the particle's center and falls off to zero at the edge of the smoothing radius.
  • Initial Function: A simple linear falloff (distance from radius).
  • Improved Functions: Cubed or squared functions are used for smoother transitions.
  • Volume Normalization: To ensure density calculations are consistent regardless of the smoothing radius, the smoothing function is normalized so its total volume is one. This is achieved by dividing the function's output by its calculated volume (e.g., $\pi r^h / 4$).

4. Density Calculation and Incompressibility

• Density Estimation: The density at any point is calculated by summing the mass of nearby particles, weighted by their influence (determined by the smoothing function).
• Incompressibility Goal: For liquids, density should be uniform. The simulation aims to correct density differences by moving particles from high-density to low-density areas.

5. Approximating Properties and Gradient Calculation

• Abstract Example: A made-up function is used to demonstrate how to approximate its values using particles and the smoothing technique.
• Property Approximation: A property at a point is approximated by summing particle properties, weighted by their mass, density, and smoothing influence.
• Core SPH Equation for Property Calculation: $A(x) = \sum_i \frac{m_i}{\rho_i} A_i W(|x - x_i|)$ Where:
  • $A(x)$ is the property at position $x$.
  • $m_i$ is the mass of particle $i$.
  • $\rho_i$ is the density of particle $i$.
  • $A_i$ is the property stored in particle $i$.
  • $W$ is the smoothing function.
  • $|x - x_i|$ is the distance between position $x$ and particle $i$.
• Gradient Calculation: To determine the direction of rapid change (essential for correcting density), the gradient is calculated.
  • Naive Approach: Finite differences using the calculate_property function (slow).
  • Optimized Approach: Utilizes the slope of the smoothing function and the direction towards the particle. This significantly speeds up computation.
  • Further Optimization: Pre-calculating and caching density values to avoid repeated computations within the gradient function.

6. Pressure Forces and Simulation Dynamics

• Target Density: A desired uniform density for the fluid.
• Pressure Calculation: A function converts density differences from the target density into a pressure value. This is noted as being more gas-like than liquid-like but is simple to implement.
• Pressure Force: Particles are moved along the pressure gradient to equalize density.
  • Force Calculation: Based on the pressure gradient.
  • Acceleration: Calculated as Force / Density (using density as mass per volume).
• Addressing Issues:
  • Division by Zero: Handled by ensuring particles don't calculate forces with themselves.
  • Incorrect Direction: A minus sign was added to the force calculation to ensure particles move towards lower density.
  • Shallow Smoothing Function Slope: Using a "spikier" smoothing function for pressure forces helps prevent particles from getting too close, improving stability.
  • Newton's Third Law: Implementing "shared pressure" (average pressure at two particles) ensures forces are equal and opposite between particles.

7. Performance Optimization: Spatial Grids

• The Bottleneck: Looping over all particles to calculate densities and forces is computationally expensive, especially when most particles are outside the smoothing radius.
• Solution: Divide space into a grid where cell size equals the smoothing radius. This limits neighbor searches to a 3x3 grid of cells around a particle.
• GPU-Friendly Approach:
  • A single array (Spatial Lookup) stores particle indices and their calculated cell keys.
  • Particles are sorted by their cell keys to group neighbors together.
  • A second array (start indices) stores the starting index for each cell key in the Spatial Lookup array.
  • This allows efficient retrieval of particles within a specific cell and its neighbors.
• Performance Improvement: This optimization dramatically increases frame rates (from ~5 FPS to ~120 FPS).

8. Advanced Techniques and Refinements

• Predictive Simulation (IGR): Predicting particle positions based on current velocity to better handle discrete time steps and upcoming situations. This significantly improved simulation behavior.
• Viscosity: Adding a term to simulate friction between fluid layers, smoothing out chaotic velocities and preventing overly energetic particles. This is inspired by the Navier-Stokes equations.
• Surface Tension Approximation: Initially, negative pressure was used to simulate a crude surface tension effect, preventing particles from splitting into small droplets. This was later replaced by a dedicated repulsive force for close particles.
• Stiffness/Jiggle: The "jelly-like" behavior is addressed by increasing the pressure multiplier (stiffness constant), but this requires more simulation steps per frame.
• Compute Shaders: The entire simulation is ported to a compute shader for parallel processing on the GPU, achieving high frame rates (e.g., 500 FPS).
  • Parallel Sorting: Implementing a bitonic merge sort algorithm for efficient sorting of particles by cell keys on the GPU.
• 3D Simulation: Extending the simulation to three dimensions by updating data structures (e.g., float3 instead of float2) and adjusting smoothing function parameters.
• Collision Handling: Transforming points and velocities to the local coordinate system of a moving/rotating bounding box for collision detection.

9. Future Improvements and Wishlist

• Stability and Performance: Further optimization for handling more particles.
• Parameter Finickiness: Making simulation parameters easier to tune for desirable results.
• Boundary Behavior: Improving how particles behave at the edges of the fluid.
• Object Interaction: Simulating how objects (like boats or ducks) interact with the fluid.
• Rendering: Developing more realistic fluid rendering techniques (e.g., ray marching).