I recently went this route. I didn’t want to set up or use Unity so I wrote my own 2D fluid simulator based on some of the same papers using Metal compute shaders (though I’d love to try again using webgpu). Sebastian’s video is great and the implementation is good. But this was a great (and fun) opportunity to look for ways to improve on it.

For starters, the way he’s doing the spatial lookup has poor cache performance, each neighbor lookup is another scattered read. Instead of rearranging an array of indices when doing the sort, just rearrange the particle values themselves. That way you're doing sequential reads for each grid cell you look for neighbors in, instead of a series of scattered reads. The performance improvement I got was about 2x, which was pretty impressive for such a simple change.

The sorting algorithm used isn’t the fastest, counting sort had much better performance for me and was simpler for me to conceptualize. It involves doing a prefix sum though, which is easy to do sequentially on the CPU but more of a challenge if you want to try keeping it on the GPU. "Fast Fixed-Radius Nearest Neighbors: Interactive Million-Particle Fluids", by Hoetzlein et al [0].

Or, if you want to keep using bitonic sort, you can take advantage of threadgroup memory to act as a sort of workspace during bitonic merge steps that are working on small enough chunks of memory. The threadgroup memory is located on the GPU die, so it has better read/write performance.

I ended up converting his pure SPH implementation to use PBF ("Position Based Fluids", Macklin et al, [1]), which is still SPH-based but maintains constant density using a density constraint solver instead of a pressure force. It seems to squeeze more stability out of each “iteration” (for SPH that’s breaking up a single frame into multiple substeps, but with PBF you can also run more iterations of the constraint solver). It’s also a whole lot less “bouncy”. One note: I had to multiply the position updates by a stiffness factor (about 0.1 in my case) to get stability, the paper doesn’t talk about this so maybe I’m doing something wrong.

The PBF paper talks about doing vorticity confinement. It’s implemented exactly as stated in the paper but I struggled for a bit to realize I could still do this in 2D. You just have to recognize that while the first cross product produces the signed magnitude of a vector pointing out of the screen, the second cross product will produce a 2D vector in the same plane as the screen. So there’s no funny business in 2D like I had originally thought. Though, you can skip vorticity confinement, the changes aren't very significant.

There’s a better (maybe a bit more expensive) method of doing surface tension/avoiding particle clustering. It behaves a lot more like fluids in real life do and avoids the “tendril-y” behavior he mentions in the video. "Versatile surface tension and adhesion for SPH fluids" by Akinci et al [2].

One of the comments on Sebastian's video mentions that doing density kernel corrections using Shepard interpolation should improve the fluid surface. I searched and found this method in a bunch of papers, including "Consistent Shepard Interpolation for SPH-Based Fluid Animation" by Reinhardt et al, [3] (I never implemented the full solution that paper proposes, though). There's kernel corrections, and then there's kernel gradient corrections, which I never got working. With the kernel corrections alone, the surface of the fluid seems to "bunch up" less when it moves, and it was pretty simple to implement. Otherwise, the surface looks a bit like a slinky or crinkling paper with particles being pushed out from the surface boundary.

I found [0] and [1] on my own but I found [2] through a thesis, "Real-time Interactive Simulation of Diverse Particle-Based Fluids" by Niall Tessier-Lavigne [4]. I also use the 2nd order integration step formula from that paper. It has some other excellent ideas that are worth trying.

Many years ago I used a paper (that is in fact one referenced by Sebastian’s video) and some C sample code I found to write an SPH simulator in OpenCL. I had been wanting to write one again but this time get a real understanding of the underlying mathematics now that I have some more tools under my belt. I owe it to Sebastian that I finally started on my implementation and I understand SPH a lot more now.

[0]: https://on-demand.gputechconf.com/gtc/2014/presentations/S41...

[1]: https://mmacklin.com/pbf_sig_preprint.pdf

[2]: https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&d...

[3]: https://www.hdm-stuttgart.de/hochschule/forschung/forschungs...

[4]: https://project-archive.inf.ed.ac.uk/ug4/20181074/ug4_proj.p...