SPH Fluid simulation

[juanchuletas]

Summary

Implement a GPU-accelerated SPH (Smoothed Particle Hydrodynamics) fluid simulation system in FunGT with screen-space fluid rendering. Target: 50K to 100K particles at interactive frame rates on Intel Arc via SYCL. Demo scene: ball dropping into a bowl of fluid.

Description

FunGT currently has a GPU rigid body physics pipeline with spatial hashing (uniform grid, radix sort, cell boundary detection). The SPH fluid simulation system will reuse this spatial infrastructure while introducing new simulation kernels for fluid dynamics and a new rendering pipeline for liquid surfaces.

Architecture: Separate SPHKernel class

The SPH system will be implemented as a new SPHKernel class, separate from the existing PhysicsKernel. Reasons:

SPH has fundamentally different per-particle data from rigid bodies. Rigid bodies carry orientation quaternions, inertia tensors, and angular velocity. Fluid particles carry density, pressure, smoothing length, and color field values. Mixing both into PhysicsKernel would bloat the DeviceData struct with fields that are never used by the other system.

The simulation pipeline is also structurally different. Rigid body physics runs: detect contacts, solve constraints, integrate. SPH runs: compute density, compute pressure, compute pressure forces, compute viscosity forces, integrate. There is no constraint solver in SPH. The forces are computed directly from the Navier-Stokes equations approximated through smoothing kernels.

The shared infrastructure (uniform grid, radix sort) can be reused by composition. SPHKernel holds a reference to the grid utilities or instantiates its own with fluid-appropriate cell sizes (cell size = smoothing radius for SPH vs cell size = 2*sphere_radius for rigid bodies).

If rigid body and fluid coupling is needed later (ball splashing into fluid), the two kernels communicate through a shared position buffer. Fluid particles detect rigid body surfaces and apply boundary forces. This is cleaner with two separate systems than with one monolithic class trying to handle both.

Rendering: Screen-Space Fluid Rendering

Real-time fluid engines do not render triangle meshes per frame. The standard pipeline (NVIDIA GDC 2010, Wicked Engine, PhysX demos) is:

1. Render particles as point sprites (one GL_POINT per particle, vertex shader computes screen-space size, fragment shader writes sphere depth and discards pixels outside the circle)
2. Output a depth texture and a thickness texture
3. Smooth the depth texture with a bilateral Gaussian blur (preserves edges while smoothing the bumpy particle surface)
4. Reconstruct surface normals from the smoothed depth texture using screen-space derivatives
5. Shade the surface with refraction, reflection, Fresnel, and thickness-based absorption in a fullscreen pass
6. Write depth to merge with the scene (fluid occludes and is occluded by solid objects)

This produces smooth, liquid-looking surfaces from raw particles without mesh reconstruction. The vertex count is exactly N particles regardless of particle count. The cost scales with screen resolution for the smoothing and shading passes, not with particle count.

Marching cubes mesh reconstruction is reserved for offline/path-traced rendering. If path-traced fluid rendering is desired later (for the FunGT path tracer), a marching cubes kernel can extract an isosurface from the SPH density field and feed it into the existing BVH-based path tracer. This is a separate task.

Particle Count: 50K to 100K

Real-time GPU SPH implementations target 10K to 1M particles depending on hardware. On Intel Arc via SYCL with OpenCL backend, 50K to 100K is a realistic target for interactive frame rates. The neighbor search using the existing radix sort and uniform grid infrastructure scales linearly with particle count. The rendering cost is dominated by the screen-space blur passes, which are resolution-dependent and constant regardless of particle count.

Requirements

Phase 1: SPH Simulation Kernel

1. New SPHKernel class with its own SPHDeviceData struct containing per-particle SoA arrays: position (x, y, z), velocity (vx, vy, vz), force (fx, fy, fz), density, pressure
2. Host-side staging (SPHHostData) for batch particle initialization, mirroring the rigid body HostData pattern
3. Uniform grid with cell size equal to the SPH smoothing radius (separate from rigid body grid)
4. computeDensity kernel: for each particle, iterate 27 neighbor cells, accumulate density contributions using the Poly6 smoothing kernel: W_poly6(r, h) = (315 / 64*pi*h^9) * (h^2 - r^2)^3
5. computePressure kernel: convert density to pressure using the equation of state: P = k * (rho - rho_rest) where k is the gas stiffness constant and rho_rest is the rest density
6. computeForces kernel: for each particle, iterate neighbors and accumulate pressure force (Spiky kernel gradient), viscosity force (Laplacian of viscosity kernel), and gravity
7. integrate kernel: semi-implicit Euler integration of velocity and position
8. Boundary handling: simple container walls (box boundaries) with penalty forces to keep fluid contained
9. Configurable simulation parameters: smoothing radius h, rest density rho_0, gas stiffness k, viscosity coefficient mu, gravity, timestep

Phase 2: Screen-Space Fluid Rendering

1. Point sprite rendering pass: draw all particles as GL_POINTS, vertex shader computes screen-space radius from world-space particle size and camera distance, fragment shader computes sphere depth and eye-space normal, outputs to depth FBO and thickness FBO
2. Depth smoothing pass: bilateral Gaussian blur on the depth texture (configurable kernel size and sigma)
3. Normal reconstruction pass: compute eye-space normals from the smoothed depth using screen-space partial derivatives of depth
4. Surface shading pass: fullscreen quad shader that computes Fresnel reflectance, refraction using the smoothed depth and normals, thickness-based color absorption (Beer's law), blends with scene background
5. Depth compositing: write fluid depth to the main depth buffer so fluid correctly occludes and is occluded by solid scene geometry

Phase 3: Demo Scene

1. Bowl or box container with visible walls
2. Fluid volume initialized inside the container (block of particles in a grid pattern)
3. Rigid body (sphere) dropped into the fluid from above
4. Fluid-rigid coupling: fluid particles detect the rigid body surface and apply boundary penalty forces; rigid body receives reaction forces from fluid (two-way coupling)

Technical Considerations

SPH Kernel Architecture

The SPHKernel class follows the same patterns as PhysicsKernel: SoA device data, staging buffers for initialization, bulk sendToDevice, SYCL kernel dispatch with DeviceData struct capture. The smoothing kernels (Poly6, Spiky, Viscosity Laplacian) are implemented as inline device functions in a header, similar to gpu_collision_detection.hpp.

The neighbor search reuses the same algorithmic pattern as detectDynamicVsDynamic: compute cell hashes with world offset, radix sort, find cell boundaries, iterate 27 neighbors. The cell size must equal the smoothing radius h so that all particles within the kernel support are guaranteed to be in adjacent cells.

Screen-Space Rendering Integration

The screen-space fluid rendering pipeline requires FBO management (depth texture, thickness texture, smoothed depth texture) and fullscreen quad passes. This integrates with the existing viewport/layer system in FunGT. The fluid rendering is a new layer or a new render pass within the existing ViewPort infrastructure.

The point sprite SSBO for particle positions uses the same SYCL-OpenCL-OpenGL interop pattern as the rigid body buildMatrices, with a persistent cl_mem object cached at init.

SPH Parameters for Water-like Fluid

Typical starting values for water-like behavior at interactive particle counts:
Smoothing radius h: 0.1 to 0.2 (depends on particle spacing)
Rest density rho_0: 1000.0 (kg/m^3 for water)
Gas stiffness k: 1000.0 to 3000.0 (higher = more incompressible, more expensive)
Viscosity mu: 0.1 to 1.0 (higher = thicker fluid like honey)
Particle mass: computed from rho_0 * volume_per_particle
Timestep: 0.001 to 0.005 (SPH requires small timesteps for stability)

Fluid-Rigid Coupling

The simplest approach is boundary particles: represent the rigid body surface as a set of fixed particles that contribute to SPH density computation but are not integrated. When a fluid particle enters the rigid body's surface, boundary particles push it out through the SPH pressure force naturally. The rigid body receives the accumulated reaction force from all fluid particles within its boundary (Newton's third law).

Acceptance Criteria

1. SPHKernel simulates 50K fluid particles at 30+ FPS on Intel Arc via SYCL OpenCL backend
2. Fluid behaves visually like water: splashes on impact, settles under gravity, fills container shape
3. Screen-space rendering produces smooth liquid surface without visible individual particles
4. Fluid correctly interacts with container boundaries (no leaking)
5. Demo scene with rigid body dropping into fluid runs at interactive frame rates
6. All SPH parameters are configurable at runtime through the existing ImGui interface
7. Point sprite SSBO uses persistent CL interop (no per-frame allocation)