July 08, 2017

A complement to my previous writeup on particle-based fluid simulation, this tutorial is meant to give a brief introduction to implementing smoothed particle hydrodynamics solvers in C++. Though I have chosen a relatively low-level compiled language for speed, these concepts should be extensible to implementations in any language, such as Python or Javascript, which might have more built-in facilities at the cost of computational speed.

A basic (unoptimized, feature-sparse) implementation of Müller's seminal SPH paper in 2D is provided on GitHub to accompany this tutorial. It is available here.

While the aforementioned writeup covers the mathematics of SPH quite well, it is worth noting once more the equations we will be directly concerned with; that is, the formulation of the incompressible Navier-Stokes equations using SPH

$\rho_i = \rho(\mathbf{r}_i)=\sum_j m_j\frac{\rho_j}{\rho_j}W(\vert\mathbf{r}_i-\mathbf{r}_j\vert,h) = \sum_j m_jW(\vert\mathbf{r}_i-\mathbf{r}_j\vert,h)\\ \mathbf{F}^{\text{pressure}}_i=-\nabla p(\mathbf{r}_i)=-\sum_jm_im_j(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2})\nabla W(\vert\mathbf{r}_i-\mathbf{r}_j\vert,h)\\ \mathbf{F}^{\text{viscosity}}_i=\eta\nabla^2\mathbf{u}(\mathbf{r}_i)=\eta\sum_jm_j\frac{\vert\mathbf{u}_j-\mathbf{u}_i\vert}{\mathbf{\rho}_j}\nabla^2W(\vert\mathbf{r}_i-\mathbf{r}_j\vert,h)$In a word, the density $\rho$ of a particle $i$ can be calculated as a weighted average over the surrounding particles; fluid forces, such as pressure and viscosity, can be modeled again summing contributions from neighboring particles, weighted by density and distance.

Keen readers will find one aspect of fluid simulation notably absent from our review discussion: surface tension. While a basic surface tension model is discussed by Müller in his original paper, code simplicity is greatly improved by focusing on modeling directly using SPH formulations.

Further, those with possible prior experience reading SPH code will note a lack of common optimizations, including, but not limited to, a griding system, parallelization, and solver term precomputation. Again, these are deliberate consequences of maintaining clean, readable, and minimal code to best explain the underlying mathematical optimizations. Explorations or tutorials regarding the optimization of SPH are planned for the future. Both griding and precomputation will be briefly touched upon in later sections of this writeup.

At the most basic level, our SPH demo must produce a reasonably realistic real-time simulation of a fluid. To this end, our code will have three principal components: a rendering system, an update system (the SPH solver), and user input handling.

We will initialize our fluid in a dam-break flow configuration, subject to boundary conditions at the edge of our rendering area, and allow the user to add additional drops of fluid.

Our rendering system will have the following methods:

```
void InitGL(void)
void Render(void)
```

Our solver:

```
void InitSPH(void)
void Integrate(void)
void ComputeDensityPressure(void)
void ComputeForces(void)
void Update(void)
```

And our interactivity handling:

```
void Keyboard(unsigned char c)
int main(int argc, char** argv)
```

Each component will be broken down in turn, with further discussion of the implementation details for each core method.

OpenGL is the gold standard for desktop graphics programming. We make use of a *very* basic OpenGL renderer built upon GLUT, a utilities toolkit which provides a cross-platform windowing system API.

The following initialization method is called by `main`

, before the rendering loop is entered. This method takes care of a few basic tasks: it sets our window background color to a nice grey, enables point smoothing (anti-aliasing), sets the point size to a multiple of our particle smoothing length, and tells the renderer we are using a view projection matrix (discussed below).

```
void InitGL(void)
{
glClearColor(0.9f,0.9f,0.9f,1);
glEnable(GL_POINT_SMOOTH);
glPointSize(H/2.f);
glMatrixMode(GL_PROJECTION);
}
```

Our render function is slightly more complicated. It is called by the windowing system each time we can redraw the current frame. That is, the effect of a real-time simulation comes from updating the simulation state, then rendering the particles in a new position, as many times per second as possible.

Breaking it down, `glClear`

gives us a blank image buffer state, upon which we can draw the contents of our new frame. `glLoadIdentity`

and `glOrtho`

set the coordinate system for drawing to be an orthographic projection in the screen space, mapped to the known constant width and height of the simulation viewing area. `glColor4f`

sets the color in RGBA for all subsequent draw calls (a nice blue for our fluid particles). Lastly, we loop through our global state array of all fluid particles, placing a vertex at the given x and y coordinate. `glBegin`

and `glEnd`

are used to signal the type of drawing we want to be applied to these vertices. We elect to place points using the flag `GL_POINTS`

, though we could just as well be placing and connecting vertices with lines (`GL_LINES`

), for example.

```
void Render(void)
{
glClear(GL_COLOR_BUFFER_BIT);
glLoadIdentity();
glOrtho(0, VIEW_WIDTH, 0, VIEW_HEIGHT, 0, 1);
glColor4f(0.2f, 0.6f, 1.0f, 1);
glBegin(GL_POINTS);
for(auto &p : particles)
glVertex2f(p.x(0), p.x(1));
glEnd();
glutSwapBuffers();
}
```

While we will not further discuss the usage of OpenGL, this type of rendering pattern is very common in simple orthographic view applications, and could easily be extended to myriad applications.

As a reminder, our solver will have the following methods:

```
void InitSPH(void)
void Integrate(void)
void ComputeDensityPressure(void)
void ComputeForces(void)
void Update(void)
```

The structure can best be seen by examining the `Update`

method first, called as frequently as possible per second between frame rendering to step our simulation state forward in time. Predictably, we follow the SPH framework of first computing densities, which we then use to compute pressure at particle locations. Density and pressure together allow us to calculate the forces acting on each particle, which are in turn used by `Integrate`

to calculate acceleration and update final particle positions for rendering.

```
void Update(void)
{
ComputeDensityPressure();
ComputeForces();
Integrate();
glutPostRedisplay();
}
```

The following data structure will track relevant quantities for all of the particles in our system: position `x`

, velocity `v`

, force `f`

, density `rho`

, and pressure `p`

.

```
struct Particle {
Particle(float _x, float _y) : x(_x,_y), v(0.f,0.f), f(0.f,0.f), rho(0), p(0.f) {}
Vector2d x, v, f;
float rho, p;
};
```

We use the fantastic `Vector2d`

provided by the Eigen library to simply our modest use of linear algebra, such as Euclidean distances.

In order to create a visually interesting simulation, we must initialize our particle system in a state that evolves interestingly in forward time. For this reason, we prefer to replicate a common setup in fluid dynamics, the so-called 'dam break', a column of fluid mass splashing down into a container of fixed size.

We achieve a 'dam break' configuration with a very simple loop, placing particles with a spacing corresponding to the smoothing length (kernel radius) of the SPH solver. The size of the dam scales as a proportion of the chosen simulation extents. A small jitter is applied to give more degenerate scattering upon contact with the floor boundary.

```
void InitSPH(void)
{
cout << "initializing dam break with " << DAM_PARTICLES << " particles" << endl;
for(float y = EPS; y < VIEW_HEIGHT-EPS*2.f; y += H)
for(float x = VIEW_WIDTH/4; x <= VIEW_WIDTH/2; x += H)
if(particles.size() < DAM_PARTICLES)
{
float jitter = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
particles.push_back(Particle(x+jitter,y));
}
}
```

The first step of our SPH solver is computing densities and pressures of the fluid at new particle positions. For all particles $i$, we loop over other particles $j$, summing their density contributions, $\rho$, subject to the application of Müller's *Poly6* kernel function. Particles beyond the smoothing radius of support, $H$, do not contribute to the summed density. The ideal gas law, $P=k_p(\rho-\rho_0)$, is chosen as an equation of state, relating pressure $p$ to the summed density.

```
void ComputeDensityPressure(void)
{
for(auto &pi : particles)
{
pi.rho = 0.f;
for(auto &pj : particles)
{
Vector2d rij = pj.x - pi.x;
float r2 = rij.squaredNorm();
if(r2 < HSQ)
{
// this computation is symmetric
pi.rho += MASS*POLY6*pow(HSQ-r2, 3.f);
}
}
pi.p = GAS_CONST*(pi.rho - REST_DENS);
}
}
```

Similarly, for all particles $i$, we loop over other particles $j$, summing their contributions to pressure and viscosity forces. Note that a simple optimization is used to precompute the constant scalar component of the gradient of the *Spiky* kernel and the Laplacian of the *Viscosity* kernel. The force of gravity is calculated from a gravitational acceleration constant `G`

and $\mathbf{F}=m\mathbf{a}$.

```
void ComputeForces(void)
{
for(auto &pi : particles)
{
Vector2d fpress(0.f, 0.f);
Vector2d fvisc(0.f, 0.f);
for(auto &pj : particles)
{
if(&pi == &pj)
continue;
Vector2d rij = pj.x - pi.x;
float r = rij.norm();
if(r < H)
{
// compute pressure force contribution
fpress += -rij.normalized()*MASS*(pi.p + pj.p)/(2.f * pj.rho) * SPIKY_GRAD*pow(H-r,2.f);
// compute viscosity force contribution
fvisc += VISC*MASS*(pj.v - pi.v)/pj.rho * VISC_LAP*(H-r);
}
}
Vector2d fgrav = G * pi.rho;
pi.f = fpress + fvisc + fgrav;
}
}
```

Whilst the equations of fluid motion and the numerical approximation approach might be new, hopefully our chosen numerical integration scheme will not be.

As SPH gives us an approximation of *density* and the *force* acting on each particle for each timestep, we can use a standard numerical integration scheme to solve for velocity, allowing us to update particle positions in forward time. At the risk of an overly verbose explanation, with known quantities of *density* and *force*, we obtain particle acceleration

We then trivially apply Euler's method, a first-order forward numerical integration scheme

$\begin{aligned} \mathbf{v}_{t+1} &= \mathbf{v}_{t} + \Delta t\cdot \mathbf{a}_{t+1} \\ &= \mathbf{v}_{t} + \Delta t\cdot \frac{\mathbf{F}_{\text{total}}}{\rho} \\ \mathbf{x}_{t+1} &= \mathbf{x}_{t} + \Delta t\cdot \mathbf{v}_{t+1} \end{aligned}$giving our final update equation for particle positions. Our `Integrate`

method applies the same to all particles:

```
void Integrate(void)
{
for(auto &p : particles)
{
// forward Euler integration
p.v += DT*p.f/p.rho;
p.x += DT*p.v;
// enforce boundary conditions
if(p.x(0)-EPS < 0.0f)
{
p.v(0) *= BOUND_DAMPING;
p.x(0) = EPS;
}
if(p.x(0)+EPS > VIEW_WIDTH)
{
p.v(0) *= BOUND_DAMPING;
p.x(0) = VIEW_WIDTH-EPS;
}
if(p.x(1)-EPS < 0.0f)
{
p.v(1) *= BOUND_DAMPING;
p.x(1) = EPS;
}
if(p.x(1)+EPS > VIEW_HEIGHT)
{
p.v(1) *= BOUND_DAMPING;
p.x(1) = VIEW_HEIGHT-EPS;
}
}
}
```

We enforce the simulation bounds by means of a simple set of boundary conditions: if a particle's position is greater than the simulation bounds (defined as `VIEW_WIDTH`

by `VIEW_HEIGHT`

), we dampen its velocity by some factor and move it back within bounds. Note that `BOUND_DAMPING`

$<0$ such that the applied velocity impulse reflects particles away from the boundary.

This unoptimized simulation loops over all neighboring particles for each particle of interest in the density and force steps; it has a solver complexity of $O(n^2)$. Thus, scaling to large numbers of particles is untenable, even with parallelization on modern CPUs. To solve this problem, many SPH implementations exploit the fact that particles beyond the kernel radius of support, $H$, do not contribute to calculations. With this in mind, a spatial grid can be constructed with cell size $2H$, meaning that any neighboring particle must exist within the directly adjacent 8 cells. Pre-sorting particles into a spatial grid for fast neighborhood iteration allows for large speedups versus unoptimized nested loops.

Due to the use of the ideal gas equation one can see that the pressure acting on a particle can quickly induce a large force with relatively little divergence from the rest density $\rho_0$, if the pressure constant $k_p$ is sufficiently large.

Further, without going into too much detail about numerical stability, one can understand fundamentally that Euler integration could exacerbate such impulses—with a large timestep, the changes in velocity and position in one frame could be well beyond the reasonable extents of the simulation, leading to equally large corrective forces the next frame, and so on.

For these reasons, our simulation is *highly* sensitive to initial parameterization, especially the gas constant, enforcing stiffness/incompressibility, and the rest density, affecting the magnitude of pressure forces. The timestep size, greatly affects the overall stability of the numerical integration and, in some rough sense, the likelihood of devolution into a chaotic mess. As an exercise to the reader, I encourage exploring the sensitivity to parameters by trying different values for these constants. I also encourage implementing an implicit integrator, such as the leapfrog method, which is stable for oscillatory systems. Forward Euler is generally not stable for these kinds of systems.

As a final note, the constant terms of the kernel functions (and their derivatives) that are used in our SPH formulation are precomputed for speed. These follow directly from analytically working out the gradient and laplacian.

```
// solver parameters
const static Vector2d G(0.f, 12000*-9.8f); // external (gravitational) forces
const static float REST_DENS = 1000.f; // rest density
const static float GAS_CONST = 2000.f; // const for equation of state
const static float H = 16.f; // kernel radius
const static float HSQ = H*H; // radius^2 for optimization
const static float MASS = 65.f; // assume all particles have the same mass
const static float VISC = 250.f; // viscosity constant
const static float DT = 0.0008f; // integration timestep
// smoothing kernels defined in Müller and their gradients
const static float POLY6 = 315.f/(65.f*M_PI*pow(H, 9.f));
const static float SPIKY_GRAD = -45.f/(M_PI*pow(H, 6.f));
const static float VISC_LAP = 45.f/(M_PI*pow(H, 6.f));
// simulation parameters
const static float EPS = H; // boundary epsilon
const static float BOUND_DAMPING = -0.5f;
```

As previously mentioned, we make use of GLUT, a cross-platform windowing toolkit, to facilitate user interaction. The user can take two simple actions by means of keypress: adding small blocks of fluid to drop and resetting the simulation. These actions are mapped to the `spacebar`

and `r`

keys, respectively. In the case of adding fluid, we use a technique similar to the *dam break* initialization, this time with a smaller block and closer spacing, giving a satisfying expansion and fall down into the container. A global counter is used to limit the number of particles.

```
void Keyboard(unsigned char c, __attribute__((unused)) int x, __attribute__((unused)) int y)
{
switch(c)
{
case ' ':
if(particles.size() >= MAX_PARTICLES)
std::cout << "maximum number of particles reached" << std::endl;
else
{
unsigned int placed = 0;
for(float y = VIEW_HEIGHT/1.5f-VIEW_HEIGHT/5.f; y < VIEW_HEIGHT/1.5f+VIEW_HEIGHT/5.f; y += H*0.95f)
for(float x = VIEW_WIDTH/2.f-VIEW_HEIGHT/5.f; x <= VIEW_WIDTH/2.f+VIEW_HEIGHT/5.f; x += H*0.95f)
if(placed++ < BLOCK_PARTICLES && particles.size() < MAX_PARTICLES)
particles.push_back(Particle(x,y));
}
break;
case 'r':
case 'R':
particles.clear();
InitSPH();
break;
}
}
```

Finally, we arrive at the `main`

method, which puts all of the pieces together. First, a window is initialized with the desired width, height, and title. `glutDisplayFunc`

sets our `Render`

function to be used to generate a new frame as often as possible. In the time between rendering, the `glutIdleFunc`

will call `Update`

, moving the simulation forward in time. `InitGL`

sets up our rendering environment, followed by `InitSPH`

, which populates our simulation with particles in a *dam break* configuration. Finally, entering the `glutMainLoop`

begins calling the render and update functions, displaying the window and blocking until it is closed by the user or the program is killed.

```
int main(int argc, char** argv)
{
glutInitWindowSize(WINDOW_WIDTH,WINDOW_HEIGHT);
glutInit(&argc, argv);
glutCreateWindow("Müller SPH");
glutDisplayFunc(Render);
glutIdleFunc(Update);
glutKeyboardFunc(Keyboard);
InitGL();
InitSPH();
glutMainLoop();
return 0;
}
```

Please see the final section of my previous tutorial for a discussion of further reading into more modern Smoothed Particle Hydrodynamics solvers. To summarize the limitations of this demonstration, surface tension has not been implemented, the lack of a gridding system severely impacts performance, no attempt was made to parallelize any computation, the ideal gas equation causes incompressibility to be poorly enforced, and the numerical integration scheme leads to parameter sensitivity and instability. With that said, I hope that this has been a useful gentle introduction to the mathematics and code of implementing real-time fluid simulations. As always, please reach out with any comments or corrections.

Last updated May 16, 2020