Wet Paint: A frame from a simulation involving 44,000 particles. Rendered in my own ray tracer.
Basics:
I first implemented almost everything in the original SPH fluids paper, Particle-Based Fluid Simulation for Interactive Applications (2003). I implemented the density and pressure computations, pressure force, viscosity force, gravity, a uniform grid with cells the size of the smoothing length for quick lookup of neighboring particles, and all of the kernel functions mentioned (poly6, spiky, and viscosity) as well as the gradients and Laplacians of some of those. I implemented symmetric force computations from another paper to comply with Newton's laws and cut the number of required computations in half.
Weakly Incompressible SPH:
Very low compressibility is an important property of water and other liquids. After looking through some more modern papers, I decided to implement weakly-compressible SPH from Weakly compressible SPH for free surface flows (2007). Weakly-compressible SPH uses a different computation for pressure (a Tait equation of state) and is more incompressible and realistic than the original SPH, but it also can require smaller time-steps because it effectively uses very stiff springs for the pressure. The paper also introduced a new and improved method for surface tension, which I also implemented.
Boundaries:
Boundaries are still an open problem in SPH fluid simulation. There are many available methods but none are perfect. A few popular methods are simply reflecting particle velocities at boundaries, using boundary particles, and using ghost-particles. One major problem with SPH is the erroneous low density at solid and free boundaries due to using an isotropic kernel to compute the density. This is a pain. In my simulations, this caused particles to clump at the boundaries to achieve the proper density, and that in turn caused the density of the inner particles to get too high and a gap formed between the outer layer of particles and the inner layers. This looked bad and caused unrealistic behavior. (Surface tension helps flatten the surface a little, but does not fix the floating layer or unrealistic behavior problems, and surface tension would round off the sharp corners and edges the simulation in a box without ghost particles.) I implemented ghost-particles for the boundaries of the box in my simulation. The ghost particles are temporary copies of only the particles that are near the boundaries. They are almost identical to the original particles, except with the normal component of velocity reflected (in the free-slip condition). I also have the necessary ghost particles at corners and edges, not just flat surfaces (figured this out, not based on any paper in particular). I also reflect the velocity (using a coefficient of restitution) if the particles actually leave the box. I also added a cylindrical boundary with ghost particles.
Surface Reconstruction:
Surface reconstruction is another important part of realistic-looking fluid simulation using SPH. I started by implementing metaballs in my own ray tracer using ray marching. The results looked good but rendering was pretty slow. So I decided to implement marching cubes to polygonize the iso-surface. But just naively using metaballs results in blobby, bumpy looking fluid—even flat surfaces look bumpy. So I looked into better solutions and found a new paper called Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic Kernels (2010) (which received the award for second-best paper at SCA 2010). This method creates elliptical kernels that do a much better job of smoothly and accurately representing the surface than existing techniques. It's important to note that the smoothing and anisotropic kernels are not propagated into the simulation—they are just used for rendering. I've implemented pretty much all of that paper including the diffusion smoothing step (which smooths out the fluid, and fixes the floating layer problem), the anisotropic kernels, the spatial data structure for storing the elliptical particles, and surface reconstruction using marching cubes. I found a few problems with the method as described in the paper that I had to work around, although it's possible that I just misinterpreted something.
Code:
I've implementing as much of the project as I can from scratch. For the OpenGL rendering, I'm using code that I wrote for another project that features a 3D scene graph, camera controls, and simple lighting. I used NIST's JAMA and TNT libraries for a singular value decomposition necessary to compute the anisotropic kernels. I'm currently using some marching cubes code from online, but have modified it to serve my purposes, in particular creating a continuous mesh instead of just a triangle soup, which is essential for shading and exporting to OBJ for rendering. I wrote the OBJ exporting myself, which was easy after I stitched the marching cubes results into a continuous triangle mesh. I used stb_image_write for outputting images, so I can play back animations.
Miscellaneous:
I have been initializing the particles in a uniform grid but giving them a random velocity at first to stir them, because uniform grids of particles can cause noticeable funny behavioral artifacts.
