Skip to content

walkthrough.md

Latest commit

 

History

History
208 lines (162 loc) · 8.24 KB

walkthrough.md

File metadata and controls

208 lines (162 loc) · 8.24 KB
Raw

SPH Fluid Simulator — Code Walkthrough

Project Structure

SPH_simulator/
├── CMakeLists.txt           # Build system (C++17, OpenGL, OpenMP)
├── include/
│   ├── SPHParticle.h        # Particle data struct (legacy, unused in SoA)
│   ├── SPHSystem.h          # SPH simulation engine header
│   ├── Shader.h             # OpenGL shader loader (header-only)
│   └── Camera.h             # FPS camera controller (header-only)
├── src/
│   ├── main.cpp             # Application entry, render loop, ImGui UI
│   ├── SPHSystem.cpp         # Core SPH physics + spatial grid
│   ├── glad.c               # OpenGL loader
│   └── imgui*.cpp           # Dear ImGui source files
├── Shader/
│   ├── particle.vert/.frag  # Particle rendering shaders
│   └── ground.vert/.frag    # Ground plane + box shaders
└── libs/
    └── glfw3.lib            # GLFW pre-built library

Build System — CMakeLists.txt

  • C++17 with MSVC optimizations (/O2 /fp:fast /openmp) enabled even in Debug
  • Links GLFW, OpenGL, and OpenMP
  • OpenMP parallelizes the SPH computation across all CPU cores

Core Physics — SPH Algorithm

Based on Müller et al. 2003 ("Particle-Based Fluid Simulation for Interactive Applications").

Data Layout — SPHSystem.h

Uses Structure of Arrays (SoA) for cache performance:

std::vector<glm::vec3> pos;       // positions
std::vector<glm::vec3> vel;       // velocities
std::vector<glm::vec3> force;     // accumulated forces
std::vector<float> density;       // computed density per particle
std::vector<float> pressure;      // computed pressure per particle

Key SPH parameters:

Parameter Default Purpose
H 0.16 Smoothing radius
GAS_CONST 1500 Pressure stiffness (k)
REST_DENSITY 20 (auto) Target rest density
VISCOSITY 8.0 Fluid viscosity (μ)
DT 0.002 Physics timestep
MASS auto Per-particle mass (calibrated)

Simulation Pipeline — SPHSystem.cpp

Each Update() call runs this pipeline:

flowchart LR
    A[BuildGridAndReorder] --> B[ComputeDensityPressure]
    B --> C[ComputeForces]
    C --> D[ApplyBoundaryForces]
    D --> E[Integrate]
Loading

1. BuildGridAndReorder()

  • Divides space into cells of size H
  • Counting sort: counts particles per cell → prefix sum → reorders particle data
  • Particles are physically rearranged in memory so neighbors are contiguous
  • Uses cellStart[]/cellEnd[] flat arrays for O(1) cell lookup

2. ComputeDensityPressure()

  • For each particle, searches 27 neighboring cells (3×3×3)
  • Density via Poly6 kernel: W(r,h) = (315/64πh⁹)(h²-r²)³
  • Pressure via gas state equation: p = k(ρ - ρ₀)
  • OpenMP parallelized with schedule(dynamic, 128)

3. ComputeForces()

  • Pressure force (Spiky gradient kernel): repels compressed particles
  • Viscosity force (Viscosity Laplacian kernel): smooths velocity differences
  • Gravity: F = ρ · g
  • Inner loops use scalar float math (not glm) for minimal overhead

4. ApplyBoundaryForces()

  • Spring-like repulsion within distance H of each wall
  • F = k_wall × (range - d) - damping × v

5. Integrate()

  • Semi-implicit Euler: velocity updated first, then position
  • Hard clamp with 0.3× velocity reversal as safety net

Kernel Coefficients (Precomputed)

POLY6_COEFF    = 315 / (64π h⁹)      // density
SPIKY_GRAD_COEFF = -45 / (π h⁶)      // pressure force
VISC_LAP_COEFF = 45 / (π h⁶)         // viscosity force

Auto-Calibration (in Reset)

Two-pass calibration ensures physics works at any particle count:

  1. Pass 1: Set MASS=1.0, measure kernel sum from particle arrangement
  2. Pass 2: Set MASS = 20.0 / kernelSum so REST_DENSITY = 20 always

Spatial Acceleration — Flat Grid

┌───┬───┬───┬───┐
│ 0 │ 1 │ 2 │ 3 │  ← cellStart[i] = first particle index in cell i
├───┼───┼───┼───┤     cellEnd[i]   = one past last particle index
│ 4 │ 5 │ 6 │ 7 │
├───┼───┼───┼───┤  Particles physically sorted by cell:
│ 8 │ 9 │10 │11 │  pos[0..2] = cell 0's particles
└───┴───┴───┴───┘  pos[3..5] = cell 1's particles, etc.
  • Grid dimensions: ceil(boxSize / H) + 2 per axis (margin cells)
  • Neighbor search: only 27 adjacent cells instead of all N particles
  • O(N×k) complexity where k ≈ 20-40 neighbors

Rendering

Particle Shader — particle.vert / particle.frag

  • Renders GL_POINTS with perspective-scaled gl_PointSize
  • Fragment shader creates sphere-like appearance:
    • Discards pixels outside unit circle
    • Fake normal from gl_PointCoord for diffuse + specular shading

Ground Shader — ground.vert / ground.frag

  • Standard Phong lighting (ambient + diffuse + specular)
  • Used for both ground plane and bounding box wireframe

Application Layer — main.cpp

Initialization

  1. GLFW window (1280×720, OpenGL 3.3 Core)
  2. GLAD loader
  3. ImGui context
  4. Shaders + SPHSystem + ground geometry

Render Loop

┌─ Process input (camera, scenario keys, pause) ──┐
│                                                  │
├─ ImGui frame (controls panel) ───────────────────┤
│                                                  │
├─ Physics update (adaptive budget, max 25ms) ─────┤
│   └─ while (accumulator >= DT && time < 25ms)   │
│       └─ sph.Update()                            │
│                                                  │
├─ Render: ground → box wireframe → particles ─────┤
│                                                  │
└─ Swap buffers ───────────────────────────────────┘

ImGui Controls

  • Scenario selection: keys 1/2/3 (one-shot)
  • Particle Count: integer input + Apply & Reset button
  • Physics: gas stiffness, viscosity, rest density, gravity, timestep, substeps
  • Rendering: point size, particle color
  • Boundary: wall stiffness, wall damping

Scenarios

Key Scenario Box Size Description
1 Blob Drop ±1 × 3 Cube of particles drops from y=1.5
2 Dam Break ±2 × 4 Wall of particles drops from y=1.5 in larger box
3 Double Blob ±1 × 3 Two blobs launched at each other

Performance Optimizations Summary

Optimization Speedup Description
Flat grid ~10-50× Replaces O(N²) with O(N×k), k≈30
Cell reorder ~2× Particles physically sorted for cache locality
SoA layout ~1.5× Separate arrays vs struct-of-arrays
OpenMP ~4-8× Parallel density + force + integration
Precomputed kernels ~1.3× No per-iteration powf() calls
Scalar inner loops ~1.2× Raw float ops instead of glm in hot path
Adaptive budget Smooth Caps physics to 25ms, drops time if slow

Dependencies

  • GLFW 3 — Windowing, input, OpenGL context
  • GLAD — OpenGL 3.3 function loader
  • GLM — Math library (vectors, matrices)
  • Dear ImGui — Immediate-mode UI for parameter tuning
  • OpenMP — CPU parallelization