This is an introduction to some of the core concepts of SPH fluids and the first part of a series on making a fluid FLIP/PIC simulation in Go. DIY fluids lets go.
For some reason since I can remember I’ve been interested in this area of research. If had to pick a reason it might be a sort of an aesthetic appeal. But mostly I think I associate fluids with a sort of visual tactile play, like its just interesting and I’m a damn child. So, I spent many hours, looking at the Navier-Stokes equations and trying to unlock its magic that probably any engineering student would know by heart. As it turns out while it’s great to understand Navier-Stokes (NS), it’s sometimes better to just get your hands dirty, knowing you can come back to NS as your guide. We won’t be doing any esoteric derivations here from Newtonian Mechanics. But here it is. Let your eyes feast.
This momentum equation really only has two key ingredients and its simple enough. A gradient pressure term and a 2nd order laplacian for viscosity. Plus the external forces. Simple enough for now although we will dive back in later chapters.
There are some other questions to ask. Like why Go? And this has been done before just search GitHub? Surely some university SIMD BLAS implementation exists. But it’s fun to make something new. And plus there weren’t any Go implementations so I figured this would be a great opportunity to add an API out to the open source world and get familiar with Go. This project is where I discovered part of why Go is so awesome actually. Although that’s maybe for another article.
In the world of computational fluid dynamics there are really two approaches, the Eulerian & Langrangian methods. I would say if you really want to learn the breadth of CFD, write a fluid solver with Grid methods and one with SPH methods. The goal for us is how to make a Hybrid solver. So we will need to cover both.
Eulerian computational methods use grids, sometimes called Finite Element/Volume Methods. Solving fluids with grids is straightforward, for the differential operators, we can implement them directly from their standard form by simple central differencing. Advection problems are handled with Semi-Langrangian methods that advect properties to other cells by backwards differencing.
An important factor when considering grids is the time-dimensional limits per step. Essentially the largest time integration step we can take in a grid solver cannot span more than two cells. This is also known as a Courant CFL condition. However, it turns out that this condition is significantly relaxed compared to the restrictions imposed by time steps for SPH. and when using Backwards Euler & Semi-Langrangian methods, we can largely overcome time step issues making grids a powerful platform for more real-time methods.
On this backwards euler do encounter one problem which is solving for diffusion due to viscosity, because viscous effects propagate across the spatial domain. To solve this case in Grids we use Backwards Euler methods which forms a linear system of equations (a vector of coefficients), whose current state can be approximated by solving a the equation Ax = b.
For this we have to process a valid LU Decomposition for all N cell parameters O(n³) typically, and carry out the Matrix-Vector Multiplication which is going to be O(n²) for our purposes. Also not that another key factor with grids is going to be handling pressure correction when implementing incompressibility, or specific compressibility conditions. This approach is handled much more succinctly in grids but still requires a solution through a linear solver.
All in all grids are great, they are reliable, solvable, and as long as your intentions are to slosh water around in a box, they are fairly easy to configure. Also if you wish to handle incompressibility you might naturally opt to solve the density and pressure solvers on a coarser grid. Not to mention Grids naturally translate to GPU shader routines since their structures are fixed.
Next up is the Langrangian framework, which in CFD is the Smoothed Particle Hydrodynamics (SPH) approach. Langrangian methods are all about moving from a spatial-framework to a framework which analyzes the system from an open set of coordinates (particles), and understanding the solution in terms of relationships between these particles. As it turns out this framework renders certain problems in FEM methods, such as arbitrarily shaped fields, collision handling, and boundary methods, implicit.
However, in SPH we no longer have a rigid static framework to track the fluid in space. In the grid methods we can track the density changes, advect some parameters according to the velocity vector field solutions, and voila, we have a convincing fluid. So essentially we can use projection methods which are reasonable to enforce constant density if that’s our goal. Most real time systems tend to just ignore incompressibility. SPH has a tendency to be a little bit “hairy”, we estimate and map density information based on neighboring particle configurations. It can be less forgiving, even if we choose to model a weakly compressible fluid.
The Kernel & EOS Equation
We achieve this density estimation through Spatial Kernel functions and there is plenty of literature on how to employ such a kernel which we will discuss later. But the main take-away is that we don’t just have incompressibility baked into the solution anymore, the fluid’s core descriptive property, it’s density, is an emergent property, not an inherent one. Eulerian methods offered us a framework with some good assurances, SPH threatens all of that with the strange and different route towards a solution.
For the SPH Kernel the idea is that we are taking a point particle in 3D space and interpolating it’s properties throughout a radial sphere. The sphere is a simple (d) distance parameter away from the particle. Now when we try to estimate that particles density we actually are taking a sum of weighted particle density contributions from other particles in the near neighborhood. (Which we track through a voxelized grid). The most common type of kernel, which is fairly flexible for most implementations is the cubic kernel.
These three kernels are derivations of the same cubic kernel where r is the distance between particles, h is the kernel smoothing length, and all other coefficients essentially ensure that the kernel integrates to 1. Kernel functions return a weighted value W for the input distance r. The first Kernel is the density kernel, which estimates the original function, based on particle distances. The Gradient & Laplacian (Pressure, Viscosity) kernels are 1st and 2nd order differentials of the original kernel function.
Grids take for granted SPH’s most perplexing and intractable problem, the incompressibility condition of the Navier-Stokes derivation. For Grid’s the density of the fluid is simply assumed for SPH it is calculated. We enforce this incompressibility in SPH by mapping an Equation-of-State (EOS) that maps the (estimated density/target density) ratio to a pressure value for that point in space. The pressure evaluation is exponential compared to the differential ratio, so that when a particles kernel sphere is filled with the appropriate configuration of neighbors, and it’s density matches the target density, the evaluated pressure at that particle is just a constant pressure, but if the density rises because of more particles in your neighborhood, well the mapped pressure skyrockets.
This EOS which has corollaries to Ideal Gas Laws, is really related to the molecular configuration of fluids, where fluid state is intractably linked to the atomic forces present, and you’d need a whole lot of pressure to even just alter those configurations for example from a liquid to a solid state under pressure. Theoretically if you could apply a sort of infinite pressure to water, after boiling hot ice, you’ll just end up invoking a fusion reaction. But in SPH typically if you have some unsatisfactory initialization, your SPH fluid is probably just destined to explode, inside your computer of course :)
In order to relax the particles into feasible states, it turns out these incredible (physics bending) forces must be resolved, typically by slowing the simulation down to unimaginably small time increments, running the risk of halting your entire simulation, otherwise risking these unrealistic (error) force terms tearing your fluid apart :(
And that’s the dark irony, SPH promises really so many benefits, and is really useful for simulating deeply accurate physics models, over a free-space, and as long as you have valid configurations. So that when body forces don’t exponentiate in any appreciable manner, you have a practical real-time method for simulating millions of particles in parallel. SPH is a double edged sword as they say. But the same goes for Grid solvers, our issue in general for water like fluids is that of handling the incompressibilty condition, a constraint that ties our hands when talking about CFD. In the next chapter we will solve for compressible fluids, then discuss incompressibility more, and see what schemes are available to cope.
Stay tuned for Part II…