Working on top of our webgl engine, Alex Miller and I created a fully gpu powered SPH, or smoothed particle hydrodynamics, fluid simulator. A demo (fully functional in chrome on most machines as of April 2016), can be found here.
To achieve the speed required for real-time frame-rates with the limitations imposed by WebGL, we had to use a few tricks.
To get all of our calculations running in parallel, we broke the steps of our simulation loop into multiple shader programs. Each program performs calculations on input which is then processed by the next stage
Since we had to limit access to shared data, each step passes information on by writing to textures.
Every time we step the simulation forward, the neighbors of every particle must be found again. Searching over every particle in the simulation was far too time-consuming, so we used a voxel grid to bin them for easier lookup.
A special shader processed each particle, and wrote its memory address into the R, G, B, or A channel of the voxel texture. The write-location in the texture corresponded to the voxel containing the physical location of the particle. This way to find a particle's neighbors, we merely had to check what particles were in surrounding voxels, and then fetch only the information about those particles from texture memory, solving an n^2 problem in linear time.
The main drawback of the voxel based approach is that it limited the accuracy of the simulation. Since each voxel was represented by only one pixel in texture memory, only four particles could be stored per voxel. Under normal circumstances this would was sufficient, but when the particles experienced large forces problems arose. The particles would be forced together, violating the law of incompressibility. Normally, a pressure force would counteract this compression, but if particles became close enough, they would start to overflow their voxels. When pressure forces should have been at their highest, overflow would cause some particles to be unaccounted for, limiting the the system's ability to restore physically accurate density. This results in a fluid that feels more squishy and compressible, and also does not transmit force as well, dampening the simulation. A potential solution would be to increase the number of pixels dedicated to each voxel in the storage texture.
To save time calculating surface tension, I came up with an algorithm that takes advantage of the voxel grid to create an approximate surface tension force with zero processing overhead. The end result of any surface tension calculation is the magnitude and direction of a force vector for each particle in the simulation. For my fake surface tension, I defined the magnitude of the vector as function of the number of empty voxels in the neighborhood of a particle. The direction of surface tension is derived from a normalized average of the vectors from the current particle to all of its neighbors. Although simple, this proved quite effective.
Another trick we came up with was a solution for dynamic collisions with objects of arbitrarily high polygon counts. From an overhead camera a depth map and a normal map are stored every frame. During the final step of our simulation, particle positions are compared with the depth map, and opposing forces are calculated based on the normal map and the distance a particle has penetrated into a scene object.
This technique has the advantage of easily handling dynamic objects, as well as very high polygon counts. It does not require creation of, or multiple lookups into, any sort of KD tree type data structure. It does however have one severe drawback. Since the undersides of objects cannot be captured in the collision texture maps, particles will not register collisions with them. This leads to some pretty un-physical behavior in certain scenes. The technique really shines when used with terrain, as most terrain models include very little or no overhanging geometry.