SPH

Welcome! I’m Owen, and this is the development blog for my 2020¹ Lagrangian² fluid simulation project. Here I’ll document progress updates, design decisions, and distillations of fluid simulation principles. Hopefully it’ll be a good time. :)

As of today, I’ve implemented a barebones SPH fluid simulation on the basis of Steve Rotenberg’s lecture slides, the 2014 SPH STAR report, this 2003 paper by Müller et al., and this report by Kelager. SPH, short for “smoothed particle hydrodynamics,” basically means that we compute fluid quantities as functions of nearby particles, where the functions vary smoothly with respect to distance. Note that in our Lagrangian frame of reference, each particle represents a small sample of the fluid, and has a particular radius of influence called a support radius.

class Particle {
    float mass, density, pressure;
    vec3 position, velocity, acceleration;
};

Preliminaries

In general, a physical simulation can be abstracted into the following loop:

while (not done)
  update:
    compute forces
    integrate motion
  draw:
    render current state

SPH is no exception. I follow Algorithm 1 from the STAR report, which describes what to do for a single update step. I also render the particles as spheres using OpenGL.

Navier-Stokes

According to the Navier-Stokes equations, the acceleration of a particle is given by

\[\begin{aligned} \mathbf{a} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{v} + \frac{\mathbf{f}_\text{other}}{m} \end{aligned}\]

Broadly, \(-\frac{1}{\rho}\nabla p\) is the acceleration due to pressure (\(\rho\) is density, \(p\) is pressure), \(\nu\nabla^2\mathbf{v}\) is the acceleration due to viscosity (\(\nu\) is the viscosity coefficient, \(\mathbf{v}\) is velocity), and \(\frac{\mathbf{f}_\text{other}}{m}\) is the acceleration due to other forces such as gravity (\(m\) is mass).

Thus, by Newton’s second law, the net force acting on a particle is

\[\begin{aligned} \mathbf{f} &= m\mathbf{a} \\ &= -\frac{m}{\rho}\nabla p + m\nu\nabla^2\mathbf{v} + \mathbf{f}_\text{other} \end{aligned}\]

Kernel Functions

To quote Steve Rotenberg’s slides, “a kernel function \(W\) represents the strength of a particle’s influence as a function of the distance from the particle.”

As per Müller et al., I use different kernels for density, pressure, and viscosity.

Density Kernel (Poly6)

\[\begin{aligned} W_d(\mathbf{x}_i - \mathbf{x}_j, h) &= \begin{cases} \frac{315}{64 \pi h^9} (h^2 - \|\mathbf{x}_i - \mathbf{x}_j\|^2)^3 & \text{if } \|\mathbf{x}_i - \mathbf{x}_j\| \leq h \\ 0 & \text{otherwise} \end{cases} \end{aligned}\]

for \(h\) the support radius.

Pressure Kernel Gradient

I use Debrun’s spiky kernel to compute pressure forces. Here, a spiky kernel is advantageous because its gradient does not vanish near the center (unlike a kernel with a flat, rounded top), causing higher repulsive forces between particles which are very close together.

\[\begin{aligned} \nabla W_s(\mathbf{x}_i - \mathbf{x}_j, h) &= \begin{bmatrix} \frac{\partial W(\mathbf{x}_i - \mathbf{x}_j, h)}{\partial x_{i, x}} \\ \frac{\partial W(\mathbf{x}_i - \mathbf{x}_j, h)}{\partial x_{i, y}} \\ \frac{\partial W(\mathbf{x}_i - \mathbf{x}_j, h)}{\partial x_{i, z}} \end{bmatrix} \\ &= \kappa(\mathbf{x}_i - \mathbf{x}_j) \end{aligned}\]

for

\[\begin{aligned} \kappa &= \begin{cases} -\frac{45 (h - \|\mathbf{x}_i - \mathbf{x}_j\|)^2}{\pi h^6 \|\mathbf{x}_i - \mathbf{x}_j\|} & \text{if } \|\mathbf{x}_i - \mathbf{x}_j\| \leq h \\ 0 & \text{otherwise} \end{cases} \end{aligned}\]

Viscosity Kernel Laplacian

\[\begin{aligned} \nabla^2 W_v(\mathbf{x}_i - \mathbf{x}_j, h) &= \begin{cases} \frac{45 (h - \|\mathbf{x}_i - \mathbf{x}_j\|)}{\pi h^6} & \text{if } \|\mathbf{x}_i - \mathbf{x}_j\| \leq h \\ 0 & \text{otherwise} \end{cases} \end{aligned}\]

Computing Forces

compute forces:
  compute neighbors
  compute density and pressure
  accumulate forces from pressure, viscosity, and gravity

Neighborhood search

For each of the following computations, we will need to know the neighborhood of each particle (i.e. particles inside of the support radius). At the beginning of each time step, I build a uniform grid hash table [dividing space into (support) x (support) x (support) cells] and use it to identify the neighbors of each particle in \(O(n)\) time. I then save those neighbors as a list attached to each particle and use them to compute forces.

Hash function

\[\begin{aligned} w_\text{grid} &= (x_\text{max} - x_\text{min}) / h \\ h_\text{grid} &= (y_\text{max} - y_\text{min}) / h \\ f(x, y, z) &= w_\text{grid}h_\text{grid}\frac{z - z_\text{min}}{h} + w_\text{grid}\frac{y - y_\text{min}}{h} + \frac{x - x_\text{min}}{h} \end{aligned}\]

where \(x_\text{min}\), \(y_\text{min}\), \(x_\text{max}\), and \(y_\text{max}\) are the parameters of the axis-aligned boundary box. Technically there are also some ceils, floors, mins, maxes, and fmodfs in there (to make sure everything gets put in the right cell), but this should give you the gist.

Density and pressure computation

For particle \(i\):

\[\begin{aligned} \rho_i &= \sum_j m_j W_d(\mathbf{x}_i - \mathbf{x}_j, h) \\ p_i &= k \left( \rho_i - \rho_0 \right) \end{aligned}\]

where the sum is over the neighborhood of particles (which should include particle \(i\) itself). \(k\) is the stiffness constant and \(\rho_0\) is the rest density.

Force due to pressure

This force will push particles from high-pressure regions to low-pressure regions. Note that the sum over neighbors should not include the particle itself, as this is a force exerted on each particle by other surrounding particles.

For particle \(i\):

\[\begin{aligned} \mathbf{f}_\text{pressure} &= -\frac{m_i}{\rho_i} \nabla p_i \\ &= -\frac{m_i}{\rho_i} \left[ \rho_i \sum_{j \neq i} m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W_s(\mathbf{x}_i - \mathbf{x}_j, h) \right] \\ &= -m_i \sum_{j \neq i} m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W_s(\mathbf{x}_i - \mathbf{x}_j, h) \end{aligned}\]

Force due to viscosity

Viscosity smooths out the velocity field, “dampening” the fluid.

For particle \(i\):

\[\begin{aligned} \mathbf{f}_\text{viscosity} &= m_i \nu \nabla^2 \mathbf{v}_i \\ &= \frac{m_i \nu}{\rho_i} \left( \sum_j \frac{m_j}{\rho_j}(\mathbf{v}_j - \mathbf{v}_i) \nabla^2 W_v(\mathbf{x}_i - \mathbf{x}_j, h) \right) \end{aligned}\]

Force due to gravity

\[\begin{aligned} \mathbf{f}_\text{gravity} &= m\mathbf{g} \end{aligned}\]

where \(g = \begin{bmatrix}0 & -9.8 & 0\end{bmatrix}^\top\) assuming an upward \(y\)-axis.

Net Force

\[\begin{aligned} \mathbf{f}_n &= \mathbf{f}_\text{pressure} + \mathbf{f}_\text{viscosity} + \mathbf{f}_\text{gravity} \end{aligned}\]

Integrating Motion

integrate motion:
  do integration
  do collision handling

I use simple forward Euler integration for now.

\[\begin{aligned} \mathbf{a}_n &= \mathbf{f}_n / m \\ \mathbf{v}_{n + 1} &= \mathbf{v}_n + \mathbf{a}_n \Delta t \\ \mathbf{x}_{n + 1} &= \mathbf{x}_n + \mathbf{v}_{n + 1} \Delta t \end{aligned}\]

CFL Condition

Optionally, I constrain the time step to be at most \(\frac{0.4h}{\|\mathbf{v}_\text{max}\|}\), so that no particle moves more than four-tenths of the support radius in one physics update.

Collision Handling

If a particle is outside of the boundary box after integration, I move it back to the collision point on the boundary, and add a scaling (proportional to how far the particle had traveled past the boundary) of the boundary’s normal to the particle’s velocity.

Rendering Current State

render current state:
  display with OpenGL
  export frames to files

Currently, I render the particles as spheres using OpenGL. To create a triangle representation of a sphere, I follow the first method from this page with sector and stack counts of 36 and 18 respectively. Also, since every particle shares the same geometry, I use instancing (glDrawElementsInstanced with instanced arrays³) to improve performance.

Current Results

I set the parameters as follows (reference):

particle mass:        0.02      kg
support radius:       0.0457    m
rest density:         998.29    kg/(m^3)
stiffness constant:   3.0       N*m
dynamic viscosity:    3.5       kg/(m*s)
delta t:              0.006     s
initial spacing:      0.02716   m

700-particle test

40,960 particles

Multiple drops, 1000 particles each

Larger-scale version of initial test

New shading, 4096 particles

I pass per-particle colors (based on densities) and normals to the shader as well.

Halloween shading, ~4000 particles

More Halloween shading, ~4000 particles

FPS

Here are my numbers on a 2013 MacBook Pro.

What’s Next?

I might not get around to updating this blog for a while, because this was all I needed to do for my class project and now I have paper deadlines to worry about. But at some point I’d like to experiment with surface extraction and high-quality rendering, GPU acceleration, guiding, and more interesting scene geometries with more sophisticated collision handling.