Hello all!

.

I'm pretty much at my wit's end with SPH. I've implemented incompressible Eulerian fluid simulations before, but this is proving to be a bit more challenging.

.

So I just want to make sure I have everything down so that I'm not mistaking everything horribly. From what I understand, the algorithm for just basic SPH is:

for each particle calculate density and pressure for each particle calculate gradient of pressure and apply it to the particle's velocity apply gravity for each particle advect each particle enforce boundary conditions

.

For calculating the density, I use the poly6 kernel:

\[ W( \vec{r} ) = \frac{315}{64 * \pi * h^9} (h^2 - |\vec{r}|^2)^3 \]

.

For calculating the pressure, I use the gradient of spikey kernel:

\[ W( \vec{r} ) = \frac{15}{\pi * h^6} (h - |\vec{r}|)^3 \]

\[ \nabla W( \vec{r} ) = - \frac{45 * (h - |\vec{r}|)^2}{\pi * h^6} * \frac{\vec{r}}{|\vec{r}|} \]

.

For calculating pressure I use \( P = k * (\rho - \rho_0) \) , although I've seen a bunch of references to a equation with an exponent of 7 floating around.

.

All of this checks out with what academic papers say, but in my current implementation, the particles immediately fall to the bottom of the screen and start squeezing together. I wish I had a video of it, but it happens so fast that OBS can't capture it.

.

I personally suspect that the gradient of the spikey kernel is is the issue, just based off of how the fluid clumps together in strange ways once all the particles fall to the bottom. I'm not applying any near-pressure, just to keep everything simple, but even without it the simulation shouldn't fall to the bottom in 0.5 seconds. Another thing I'm suspicious of is the k, rho_0, and h constants, because it's a little difficult to tune the simulation when it doesn't last very long

.

Does anyone have any experience with SPH who could lend some insight? I implemented this paper but it can't be implemented very easily on the GPU, which is why I'm interested in a more pure SPH approach, so insight on either one of these would be enormously helpful!

.

I'm attaching the code for anyone who could take the time to look over it. It's 1 file, a couple functions, nothing too big or too fancy fancy of data structures. All computations are O(n^2) and structures are not very optimized, just as a sanity-check.

.

Thanks in advance!