SPH Pressure Force and Viscosity Integration in GLSL
This article provides a technical breakdown of how Smoothed Particle Hydrodynamics (SPH) models fluid dynamics on the GPU, focusing specifically on the mathematical formulations of the pressure gradient and viscosity forces derived from the Navier-Stokes equations. It covers the symmetric discretization schemes required to satisfy Newton's third law, the standard kernel functions used for numerical stability, and the translation of these continuous physics equations into high-performance OpenGL Shading Language (GLSL) compute shaders for real-time integration.
Mathematical Formulation
Smoothed Particle Hydrodynamics discretizes a continuous fluid into interpolating particles. The movement of these particles is governed by the momentum conservation equation in the Navier-Stokes model:
\[\rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \mu \nabla^2 \mathbf{v} + \mathbf{f}^{\text{ext}}\]
Where \(\rho\) is the fluid density, \(\mathbf{v}\) is the velocity vector, \(p\) is scalar pressure, \(\mu\) is dynamic viscosity, and \(\mathbf{f}^{\text{ext}}\) represents external forces like gravity.
1. Pressure Force Formulation
A direct gradient evaluation of pressure produces asymmetrical forces between particle pairs, violating the conservation of momentum. To ensure anti-symmetric forces (\(\mathbf{f}_{ij} = -\mathbf{f}_{ji}\)), the SPH pressure gradient is symmetrized:
\[\mathbf{f}_i^{\text{pressure}} = -\frac{1}{\rho_i} \nabla p(\mathbf{x}_i) = -\sum_{j} m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W(\mathbf{x}_i - \mathbf{x}_j, h)\]
Alternatively, Müller's computationally efficient formulation is frequently utilized in real-time graphics:
\[\mathbf{f}_i^{\text{pressure}} = -\sum_{j} m_j \frac{p_i + p_j}{2 \rho_j} \nabla W_{\text{spiky}}(\mathbf{r}_{ij}, h)\]
Where \(\mathbf{r}_{ij} = \mathbf{x}_i - \mathbf{x}_j\), \(r = \|\mathbf{r}_{ij}\|\), and \(h\) is the smoothing radius. Pressure \(p_i\) is derived using the Tait equation of state:
\[p_i = k \left( \left( \frac{\rho_i}{\rho_0} \right)^\gamma - 1 \right) \quad \text{or linearly} \quad p_i = k (\rho_i - \rho_0)\]
Here, \(k\) controls stiffness, \(\rho_0\) is the rest density, and \(\gamma \approx 7\) for water. To avoid particle clustering under high pressure, the Müller Spiky kernel gradient is used:
\[\nabla W_{\text{spiky}}(\mathbf{r}_{ij}, h) = -\frac{45}{\pi h^6} (h - r)^2 \frac{\mathbf{r}_{ij}}{r}\]
2. Viscosity Force Formulation
Viscosity dissipates relative kinetic energy into heat. Evaluating the Laplacian of the velocity field directly can cause numerical instabilities. The standard symmetrical formulation is:
\[\mathbf{f}_i^{\text{viscosity}} = \frac{\mu}{\rho_i} \nabla^2 \mathbf{v}(\mathbf{x}_i) = \mu \sum_{j} m_j \frac{\mathbf{v}_j - \mathbf{v}_i}{\rho_j} \nabla^2 W_{\text{visc}}(\mathbf{r}_{ij}, h)\]
The viscosity kernel designed by Müller preserves a positive Laplacian everywhere inside the support radius:
\[\nabla^2 W_{\text{visc}}(\mathbf{r}_{ij}, h) = \frac{45}{\pi h^6} (h - r)\]
Temporal Integration Scheme
Once pressure and viscosity forces are computed, net particle acceleration \(\mathbf{a}_i\) is integrated using a Symplectic Euler or Velocity Verlet scheme:
\[\mathbf{a}_i = \frac{\mathbf{f}_i^{\text{pressure}}}{\rho_i} + \frac{\mathbf{f}_i^{\text{viscosity}}}{\rho_i} + \mathbf{g}\]
\[\mathbf{v}_i(t + \Delta t) = \mathbf{v}_i(t) + \mathbf{a}_i \Delta t\]
\[\mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + \mathbf{v}_i(t + \Delta t) \Delta t\]
GLSL Compute Shader Implementation
The following GLSL compute pass calculates the pressure and viscosity forces for a given particle and updates its velocity and position.
#version 450 core
layout(local_size_x = 256) in;
struct Particle {
vec4 position; // xyz: pos, w: density
vec4 velocity; // xyz: vel, w: pressure
};
layout(std430, binding = 0) readonly buffer ParticleIn {
Particle particlesIn[];
};
layout(std430, binding = 1) writeonly buffer ParticleOut {
Particle particlesOut[];
};
uniform int uNumParticles;
uniform float uH; // Smoothing radius
uniform float uMass; // Constant particle mass
uniform float uViscosity; // Dynamic viscosity coefficient (mu)
uniform float uDeltaTime;
uniform vec3 uGravity;
const float PI = 3.14159265359;
void main() {
uint i = gl_GlobalInvocationID.x;
if (i >= uNumParticles) return;
vec3 pos_i = particlesIn[i].position.xyz;
vec3 vel_i = particlesIn[i].velocity.xyz;
float rho_i = particlesIn[i].position.w;
float p_i = particlesIn[i].velocity.w;
vec3 f_pressure = vec3(0.0);
vec3 f_viscosity = vec3(0.0);
float h2 = uH * uH;
float spiky_const = -45.0 / (PI * pow(uH, 6.0));
float visc_const = 45.0 / (PI * pow(uH, 6.0));
// In production, neighbor queries should utilize spatial hashing/grids
for (int j = 0; j < uNumParticles; ++j) {
if (i == j) continue;
vec3 pos_j = particlesIn[j].position.xyz;
vec3 r_vec = pos_i - pos_j;
float r2 = dot(r_vec, r_vec);
if (r2 > 0.0 && r2 < h2) {
float r = sqrt(r2);
float h_minus_r = uH - r;
vec3 dir = r_vec / r;
vec3 vel_j = particlesIn[j].velocity.xyz;
float rho_j = particlesIn[j].position.w;
float p_j = particlesIn[j].velocity.w;
// Pressure force: Spiky Kernel Gradient
vec3 grad_W = spiky_const * h_minus_r * h_minus_r * dir;
f_pressure -= uMass * ((p_i + p_j) / (2.0 * rho_j)) * grad_W;
// Viscosity force: Viscosity Kernel Laplacian
float laplace_W = visc_const * h_minus_r;
f_viscosity += uViscosity * uMass * ((vel_j - vel_i) / rho_j) * laplace_W;
}
}
// Combine forces and integrate
vec3 acceleration = (f_pressure + f_viscosity) / rho_i + uGravity;
vec3 new_velocity = vel_i + acceleration * uDeltaTime;
vec3 new_position = pos_i + new_velocity * uDeltaTime;
// Simple boundary condition (bounding box at y = 0)
if (new_position.y < 0.0) {
new_position.y = 0.0;
new_velocity.y *= -0.5; // Restitution damping
}
particlesOut[i].position = vec4(new_position, rho_i);
particlesOut[i].velocity = vec4(new_velocity, p_i);
}