SPH Fluid Density Evaluation in GLSL

This article explains how to compute fluid density at particle positions in Smoothed Particle Hydrodynamics (SPH) using the OpenGL Shading Language (GLSL). In SPH simulations, density estimation is the foundational step required to calculate fluid pressure and internal forces. Below, you will explore the mathematical formulation of SPH density estimation, the role of spatial data structures for efficient neighbor searching, and a concrete GLSL compute shader implementation for evaluating density in real time.

The Mathematical Foundation

Smoothed Particle Hydrodynamics approximates a continuous fluid field using a finite set of moving interpolation points (particles). The fluid density \(\rho_i\) at particle \(i\) is calculated by summing the mass contributions of all adjacent particles \(j\) within a finite interaction radius \(h\), weighted by a smoothing kernel function \(W\):

\[\rho_i = \sum_{j} m_j W(\mathbf{r}_i - \mathbf{r}_j, h)\]

Where:

  • \(m_j\) is the mass of particle \(j\).
  • \(\mathbf{r}_i - \mathbf{r}_j\) is the displacement vector between particle \(i\) and neighbor \(j\).
  • \(W(\mathbf{r}, h)\) is a radially symmetric smoothing kernel that drops to zero when the distance \(\|\mathbf{r}\| \ge h\).

A standard kernel choice for density evaluation is the Poly6 kernel due to its computational efficiency and smooth falloff:

\[W_{\text{poly6}}(r, h) = \frac{315}{64 \pi h^9} (h^2 - r^2)^3 \quad \text{for } 0 \le r \le h\]

This formulation inherently includes the particle's self-contribution when \(j = i\) (\(r = 0\)).

Spatial Acceleration via Uniform Grids

A naive evaluation over all \(N\) particles scales at \(O(N^2)\), which is prohibitive for real-time applications. To parallelize this on the GPU with GLSL, space is divided into a regular 3D grid where each cell has dimensions equal to or slightly larger than the support radius \(h\).

Before running the density compute shader:

  1. Particles are assigned a 3D cell coordinate \(\mathbf{c} = \lfloor \mathbf{r} / h \rfloor\).
  2. A spatial hash or 1D index is assigned to each particle.
  3. Particle arrays are sorted by cell ID, or a linked list/index table is constructed using Shader Storage Buffer Objects (SSBOs).

This reduces the search space for any particle \(i\) to only its current cell and the 26 surrounding neighboring cells in 3D (or 8 in 2D).

GLSL Compute Shader Implementation

The following GLSL compute shader reads particle positions from an SSBO, queries adjacent cells from an acceleration structure, evaluates the Poly6 kernel, and writes the calculated density into an output buffer.

#version 430 core

layout(local_size_x = 256, local_size_y = 1, local_size_z = 1) in;

// Particle buffers
layout(std430, binding = 0) readonly buffer PositionBuffer {
    vec4 positions[]; // .xyz = position, .w = particle mass
};

layout(std430, binding = 1) writeonly buffer DensityBuffer {
    float densities[];
};

// Spatial grid lookup buffers
layout(std430, binding = 2) readonly buffer GridCellOffsets {
    uint cellStartIndices[]; // Starting index of particles in each cell
};

layout(std430, binding = 3) readonly buffer GridCellCounts {
    uint cellParticleCounts[]; // Number of particles in each cell
};

// Simulation parameters
uniform float smoothingRadius;
uniform float poly6Coeff; // 315.0 / (64.0 * PI * pow(smoothingRadius, 9))
uniform ivec3 gridDimensions;
uniform vec3 worldBoundsMin;
uniform uint totalParticles;

// Maps a 3D grid coordinate to a 1D cell index
int getCellIndex(ivec3 cell) {
    if (any(lessThan(cell, ivec3(0))) || any(greaterThanEqual(cell, gridDimensions))) {
        return -1;
    }
    return cell.x + cell.y * gridDimensions.x + cell.z * (gridDimensions.x * gridDimensions.y);
}

void main() {
    uint index = gl_GlobalInvocationID.x;
    if (index >= totalParticles) {
        return;
    }

    vec3 pos_i = positions[index].xyz;
    float h2 = smoothingRadius * smoothingRadius;
    float density = 0.0;

    // Determine current cell coordinate
    ivec3 currentCell = ivec3(floor((pos_i - worldBoundsMin) / smoothingRadius));

    // Iterate through neighboring cells (3x3x3 volume)
    for (int z = -1; z <= 1; ++z) {
        for (int y = -1; y <= 1; ++y) {
            for (int x = -1; x <= 1; ++x) {
                ivec3 neighborCell = currentCell + ivec3(x, y, z);
                int cellIdx = getCellIndex(neighborCell);

                if (cellIdx == -1) continue;

                uint start = cellStartIndices[cellIdx];
                uint count = cellParticleCounts[cellIdx];

                // Iterate through particles in the neighbor cell
                for (uint p = 0; p < count; ++p) {
                    uint neighborParticleIdx = start + p;
                    vec4 neighborData = positions[neighborParticleIdx];
                    vec3 pos_j = neighborData.xyz;
                    float mass_j = neighborData.w;

                    vec3 diff = pos_i - pos_j;
                    float dist2 = dot(diff, diff);

                    // Apply Poly6 kernel if within smoothing radius
                    if (dist2 < h2) {
                        float term = h2 - dist2;
                        density += mass_j * poly6Coeff * term * term * term;
                    }
                }
            }
        }
    }

    densities[index] = density;
}

Critical Performance Considerations

  1. Precomputed Coefficients: Compute constants such as \(\frac{315}{64 \pi h^9}\) on the host CPU and pass them as uniforms rather than recalculating them per particle inside the shader.
  2. Squared Distances: Use squared Euclidean distances (dist2 < h2) to avoid expensive square root (sqrt) operations inside neighbor iteration loops.
  3. Memory Coalescing: Memory access patterns significantly impact GLSL execution speed. Sorting particle position and property buffers by spatial cell order guarantees contiguous memory reads when threads evaluate particles in the same neighborhood.