Smoothed-particle Hydrodynamics Fluid Simulation
1. Implement basic SPH algorithm
The basic SPH algorithm is implemented Following the general overview at SPH Fluids in Computer Graphics. SPH uses small moving particles to compute various quantities and derivatives of fluid for simulation.
Setup
I use the standard smoothing kernel described in the paper. \[W_{ij}=W(\frac{||\textbf{x}_i-\textbf{x}_j||}{h})=W(q)=\frac{1}{h^d}f(q)\] where \(d=3\), \[f(q)=\frac{3}{2\pi}\begin{cases} \frac{2}{3}-q^2+\frac{1}{2}q^3 & 0\leq q < 1 \\ \frac{1}{6}(2-q)^3 & 1\leq q<2 \\ 0 & q\geq 2 \end{cases}\]
I will use SI units through this work. For external forces, I use gravitational acceleration \([0,-9.8,0]\). I set the fluid rest density to \(\rho_0=1000\) to match the density of water. To compute pressure, I use
\[p_i = \frac{k\rho_0}{\gamma}((\frac{\rho_i}{\rho_0})^\gamma - 1) \] where \(k = 7, \gamma = 7\).
At each time step, CFL condition is tested to set the minimum time step. I use \[\Delta t = 0.2\frac{h}{||\textbf{v}^{max}||}\]
I also clamp \(\Delta t\) in the range \([0.0001, 0.003]\) since extremely small time step indicates extremely high speed caused by instability, and we need to improve the solver instead of lowering the time step. \(\Delta t\) greater than \(0.003\) does not have visible computation impact on my system so I set it as the maximum.
The viscosity computation is exactly the same as in the paper. A simple forward Euler integrator is applied each simulation step. All above are implemented with CUDA.
I handle the simulation boundary by applying a force field pointing inside. Fluid particles near the boundary will be pushed inwards. Fluid that moves out of the boundary will be moved back into the boundary.
Visualization
I visualize fluid with low and high viscosity in the following pictures. I use OpenGL 4 for the rendering.
2. Better Solver: PCI
I implemented Predictive-Corrective Incompressible SPH by Solenthaler and Pajarola. Instead of computing pressure directly from density, they initialize pressure to \(0\). Then they predict velocities, positions, and densities of the next time step. Next they solve for a pressure that corrects the density error from the predicted density and rest density.
Here is a demo of this method
3. Performance improvements
I implemented all parts of this project in CUDA. One slowest part is finding neighbors. I deployed the method described in Fast Fixed-Radius Nearest Neighbors: Interactive Million-Particle Fluids by Hoetzlein at NVIDIA. They use the counting sort algorithm on GPU to speed up the uniform grid neighbor search. This method has several steps.
- Count particles in each space cell. This is done by iterating all particles and use the atomic add operations to compute the particle count for each cell.
- Getting the sorted index of each particle. We first perform the exclusive scan (prefix sum) on the grid to get the starting index for the first particle in each cell. Then we insert particle indices into the index list.
Now each particles knows its cell and each cells knows the staring index in the particle index list. We can find neighbors efficiently.
4. Marching cube surface extraction
I implemented the marching cube algorithm with CUDA. I set the marching cube cell size to half the size of particles. I reuse \(f(q)\), accumulating the value on grid points. I use a threshold of \(0.2\) to extract the surface. For normal vectors, I compute the (negative) gradient of the accumulated values on grid points. I interpolate the grid values and normal vectors to extract a smooth fluid surface. Here are short videos for the extracted surface.
5. Ray tracer
I implemented a ray tracer with NVIDIA OptiX for real time rendering. So for my SPH pipeline, after data initialization, all simulation, surface extraction, and rendering occurs on GPU, sharing the underlying data buffers. First I present a reflection only rendering without global illumination. The rendering is close to real time. The physical simulation is slower but still at interactive speed.
Next, if we sample more rays per pixel, we can render refraction with reflection. (only rendering refraction does not look good) Here is a short video showing this. I trace 9 rays per pixel. I use 1 point light. If the camera ray hits a diffuse surface, 1 shadow ray is shot and the camera ray is terminated (no secondary reflection). If the camera ray hits a reflective or refractive surface, a new ray is generated based on Fresnel equations and a maximum of 7 bounces is allowed between these surfaces. Since the refraction of shadow ray would require other techniques to achieve (e.g. bidirectional path tracing), I do not generate physically correct water shadows and caustics. Instead, I make water attenuate the passing straight through it.
6. Spray, Foam, Bubbles
I implemented the special visual particles described in Unified Spray, Foam and Bubbles for Particle-Based Fluids. Water becomes diffuse white color when air is mixed with it. Air usually mixes with water at places where water impacts happen or at wave crest with a negative curvature. Also, more air is mixed when water moves faster.
I implemented the method described in this paper to spawn new "visual particles" at wave crests and impact sites. I also set the life time of the visual particles to be quite short so the effect is more visible in my case. The new diffuse particles are rendered as small solid Lambertian tetrahedrons that do not cast shadow. Here is a videos showing this effect.
View this project on Github https://github.com/fbxiang/fluid