Implementing SPH in 2D | Lucas V. Schuermann

Implementing SPH in 2D
Published 07/08/2017

=========================

Table of Contents

Introduction
Source Code
Review
Simplifying Modifications
Implementation
Next Steps
Citation

## Introduction

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 (SPH) 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.

## Source Code

A basic (unoptimized, feature-sparse) implementation of Müller's seminal SPH paper in 2D is provided on GitHub to accompany this tutorial. Please see the code repository (C++11), or a more recent implementation in Rust.

SPH Simulation Screenshot

## Review

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.

## Simplifying Modifications

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 mathematics. 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.

## Implementation

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.

### Rendering

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.f, 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.

### Solver

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();
}

#### Design

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 simplify our modest use of linear algebra, such as Euclidean distances.

### Initialization

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.

Dam Break Simulation Screenshot

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));
            }
            else
            {
                return;
            }
        }
    }
}

### Density and Pressure

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);
    }
}

### Calculating Forces

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, 3.f);
                // compute viscosity force contribution
                fvisc += VISC * MASS * (pj.v - pi.v) / pj.rho * VISC_LAP * (H - r);
            }
        }
        Vector2d fgrav = G * MASS / pi.rho;
        pi.f = fpress + fvisc + fgrav;
    }
}

### Numerical Integration

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

\mathbf{a}_{t+1} = \mathbf{\ddot{x}} = \frac{\mathbf{F}_{\text{total}}}{\rho}

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.f)
        {
            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.f)
        {
            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;
        }
    }
}

#### Enforcing Boundary Conditions

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.

### Performance Considerations

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.

#### A Word On Constants

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, -10.f);   // external (gravitational) forces
const static float REST_DENS = 300.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 = 2.5f;        // assume all particles have the same mass
const static float VISC = 200.f;       // viscosity constant
const static float DT = 0.0007f;       // integration timestep

// smoothing kernels defined in Müller and their gradients
// adapted to 2D per "SPH Based Shallow Water Simulation" by Solenthaler et al.
const static float POLY6 = 4.f / (M_PI * pow(H, 8.f));
const static float SPIKY_GRAD = -10.f / (M_PI * pow(H, 5.f));
const static float VISC_LAP = 40.f / (M_PI * pow(H, 5.f));

// simulation parameters
const static float EPS = H; // boundary epsilon
const static float BOUND_DAMPING = -0.5f;

### Interactivity

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;
}

## Next Steps

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 griding 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.

## Citation

Cited as:

Schuermann, Lucas V. (Jul 2017). Implementing SPH in 2D. Writing. https://lucasschuermann.com/writing/implementing-sph-in-2d

@misc{schuermann2017sph2d,
  title   = {Implementing SPH in 2D},
  author  = {Schuermann, Lucas V.},
  year    = {2017},
  month   = {Jul},
  url     = "https://lucasschuermann.com/writing/implementing-sph-in-2d"
}