Smoothed Particle Hydrodynamics
I realized that the problems I was having with the tracking of the water volume were not as easy to fix as I thought. It seemed that grid-based (Eulerian) methods are very suitable for a continuous fluid, but not so good when a sharp boundary between water and air is needed.
So I turned to particle-based (Lagrange) fluids instead. The standard way to do this is the method of Smoothed Particle Hydrodynamics, or SPH. The problem of a particle method is that you know the density, velocity etcetera only at a very limited number of places: at the particles themselves. Interpolation is not as easy as with a grid-based method, so we use a smoothing kernel over the neighbouring particles instead.
Since the number of interactions between particles is O(n^2), we disregard the influence of faraway particles, so we can use a spatial data structure and keep our simulation O(n)-ish. In my case, I added a simple grid structure and only counted the influence of particles in cell itself and its eight direct neighbours.
It took about two hours to throw together a simulation based on Müller’s work from 2003. It took about five milliseconds for that simulation to blow up.
It turns out that the biggest problem with SPH is to get the parameters right. Unlike with a grid method, most don’t even have a physical interpretation! The most important parameter is without doubt the size of the smoothing kernel. Too small, and your particles won’t affect each other anymore; values will be based on just a single particle. Too large, and you will get too much smearing of values. It’s the same with the ‘spring’ that is added between particles as a result of the pressure difference: too weak, and the particles will happily bunch together at the bottom of the volume; too stiff, and the simulation will blow up.
Even the parameters that should have a physical interpretation, don’t. At least, that’s my experience. I started with values that made sense in the real world: metres, seconds, kilograms, Pascals. These values worked well in my grid-based solver. In the particle solver, it took hours of tweaking and fiddling. In a box of size 1 by 1, gravity ended up at 10000, viscosity at 0.0001, the time step at 0.000005. But then I finally had something that slightly resembled a fluid.
It does not look like a decent water surface, because I haven’t added the surface tension term yet. That’s something for tomorrow.
But all this trouble is giving me doubts about the use of the SPH method. It’s like walking a tightrope, and the slightest change in parameter value causes the simulation to blow up. Grid methods seem much more predictable in that respect; yes, they can blow up, but at least there’s theory that tells you when you’re safe.
Hybrid approaches exist, that use a grid method to compute the flow, but particles to track the surface. They seemed needlessly complex, but they might give the best of both worlds. Or the worst. I’ll have a look tomorrow.