Lattice Boltzmann Method in GLSL Compute Shaders

This article provides a comprehensive breakdown of formulating and implementing the Lattice Boltzmann Method (LBM) for fluid flow simulation using OpenGL Shading Language (GLSL) compute shaders. It covers the core mathematical formulation of the standard D2Q9 lattice and Bhatnagar-Gross-Krook (BGK) collision operator, details the mapping of particle distribution functions to GPU memory, explains the pull-streaming paradigm essential for compute execution, and demonstrates how to handle macroscopic moment reconstruction and boundary conditions natively in shader code.

1. Mathematical Formulation: The D2Q9 LBM Model

The Lattice Boltzmann Method models fluid dynamics at a mesoscopic scale by tracking discrete velocity distribution functions \(f_i(\mathbf{x}, t)\), representing the probability of finding a fluid particle at position \(\mathbf{x}\) and time \(t\) moving with lattice velocity \(\mathbf{c}_i\).

In a standard two-dimensional, nine-velocity model (D2Q9), space is discretized into a square Cartesian grid, and particles move along nine discrete directions:

\[\mathbf{c}_0 = (0, 0)\] \[\mathbf{c}_{1,2,3,4} = (\pm 1, 0), (0, \pm 1)\] \[\mathbf{c}_{5,6,7,8} = (\pm 1, \pm 1)\]

Each direction has an associated lattice weight \(w_i\):

\(w_0 = \frac{4}{9}\)
\(w_{1..4} = \frac{1}{9}\)
\(w_{5..8} = \frac{1}{36}\)

The speed of sound on this lattice is defined as \(c_s = \frac{1}{\sqrt{3}}\), making \(c_s^2 = \frac{1}{3}\).

The BGK Collision and Streaming Equations

The discrete Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) approximation splits the evolution into two continuous steps: collision and streaming.

\[f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) - \frac{1}{\tau} \left( f_i(\mathbf{x}, t) - f_i^{\text{eq}}(\mathbf{x}, t) \right)\]

Where:

\(\tau\) is the dimensionless relaxation time, directly related to the kinematic shear viscosity \(\nu\) by \(\nu = c_s^2 \left(\tau - \frac{1}{2}\right)\).
\(f_i^{\text{eq}}\) is the Maxwell-Boltzmann equilibrium distribution function, truncated to second order:

\[f_i^{\text{eq}}(\rho, \mathbf{u}) = w_i \rho \left( 1 + \frac{\mathbf{c}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{u})^2}{2 c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2 c_s^2} \right)\]

Macroscopic Quantities

The fluid macroscopic density \(\rho\) and momentum density \(\rho \mathbf{u}\) are recovered through the velocity moments of the distributions:

\[\rho = \sum_{i=0}^{8} f_i\]

\[\mathbf{u} = \frac{1}{\rho} \sum_{i=0}^{8} f_i \mathbf{c}_i\]

2. Adapting LBM for GLSL Compute Shaders

Translating LBM to GLSL requires mapping the lattice structures to GPU memory layout and restructuring memory operations to avoid race conditions.

Memory Layout: Textures vs. SSBOs

For D2Q9, nine scalar distribution values must be stored per grid cell:

2D Texture Arrays (sampler2DArray / image2DArray): Allocate an array of nine layers of size \(W \times H\) using R32F format.
Shader Storage Buffer Objects (SSBOs): Flat arrays where index calculations follow index = (y * Width + x) * 9 + i (Array of Structures) or index = i * (Width * Height) + (y * Width + x) (Structure of Arrays).

Structure of Arrays (SoA) is preferred on GPUs because it ensures coalesced global memory access. Using three RGBA32F textures (storing 12 channels total, leaving 3 unused) or a 9-layer image2DArray of format R32F provides optimal cache locality for 2D setups.

Pull vs. Push Streaming

Standard LBM defines streaming as a push operation: a cell computes post-collision states and writes them out to its neighbors (\(\mathbf{x} + \mathbf{c}_i\)). In parallel compute shaders, pushing leads to write hazards where multiple threads attempt to write to the same memory cell.

To make the algorithm embarrassingly parallel, the formulation is inverted to a pull streaming scheme:

A compute thread is dispatched for every grid coordinate \(\mathbf{x}\).
The thread identifies its neighbors along the incoming velocity directions (\(\mathbf{x} - \mathbf{c}_i\)).
It fetches the distributions from those neighbors (Streaming).
It performs the collision locally and writes the result directly to its own cell coordinate \(\mathbf{x}\).

This requires a ping-pong buffer strategy: reading from buffer \(A\) at step \(t\) and writing to buffer \(B\) for step \(t+1\).

3. GLSL Compute Shader Implementation

Below is the complete algorithmic implementation of a combined pull-stream and BGK collision step inside a GLSL compute shader.

#version 460 core

layout (local_size_x = 16, local_size_y = 16) in;

// Ping-Pong image arrays for D2Q9 distributions (Layer 0 to 8)
layout (r32f, binding = 0) uniform readonly  image2DArray u_DistributionsIn;
layout (r32f, binding = 1) uniform writeonly image2DArray u_DistributionsOut;
// Flag image: 0 = Fluid, 1 = Solid Boundary
layout (r8ui,  binding = 2) uniform readonly  uimage2D      u_Geometry;

// Simulation parameters
uniform float u_Tau; // Relaxation time

// D2Q9 Constants
const ivec2 c[9] = ivec2[](
    ivec2( 0,  0),
    ivec2( 1,  0), ivec2( 0,  1), ivec2(-1,  0), ivec2( 0, -1),
    ivec2( 1,  1), ivec2(-1,  1), ivec2(-1, -1), ivec2( 1, -1)
);

const float w[9] = float[](
    4.0 / 9.0,
    1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0,
    1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0
);

const int opposite[9] = int[](0, 3, 4, 1, 2, 7, 8, 5, 6);

void main()
{
    ivec2 coords = ivec2(gl_GlobalInvocationID.xy);
    ivec2 size = imageSize(u_Geometry);

    if (coords.x >= size.x || coords.y >= size.y) {
        return;
    }

    uint cellType = imageLoad(u_Geometry, coords).r;

    // Solid cells do not compute collisions
    if (cellType == 1u) {
        return;
    }

    float f[9];
    float rho = 0.0;
    vec2 momentum = vec2(0.0);

    // 1. Pull Streaming Step
    for (int i = 0; i < 9; ++i) {
        ivec2 sourceCoords = coords - c[i];
        
        // Check if the source location is a solid wall
        uint sourceCell = imageLoad(u_Geometry, sourceCoords).r;

        if (sourceCell == 1u) {
            // Half-way bounce-back: reflect from boundary
            f[i] = imageLoad(u_DistributionsIn, ivec3(coords, opposite[i])).r;
        } else {
            // Stream distribution from neighboring fluid cell
            f[i] = imageLoad(u_DistributionsIn, ivec3(sourceCoords, i)).r;
        }

        // Compute zero-th and first moments
        rho += f[i];
        momentum += f[i] * vec2(c[i]);
    }

    // 2. Macroscopic velocity
    vec2 u = momentum / rho;

    // 3. Collision Step (BGK)
    float omega = 1.0 / u_Tau;
    float uSq = dot(u, u);

    for (int i = 0; i < 9; ++i) {
        float c_dot_u = dot(vec2(c[i]), u);
        
        // Maxwell-Boltzmann equilibrium
        float feq = w[i] * rho * (
            1.0 + 
            3.0 * c_dot_u + 
            4.5 * (c_dot_u * c_dot_u) - 
            1.5 * uSq
        );

        // Relax toward equilibrium and save to output buffer
        float f_collided = f[i] - omega * (f[i] - feq);
        imageStore(u_DistributionsOut, ivec3(coords, i), vec4(f_collided, 0.0, 0.0, 0.0));
    }
}

4. Boundary Handling and Numerical Stability

Bounce-Back Scheme

Rigid no-slip boundaries are handled via the standard half-way bounce-back scheme. When a streaming path crosses an obstacle node, the incoming momentum vector is flipped:

\[f_{\text{opposite}(i)}(\mathbf{x}, t + \Delta t) = f_i^*(\mathbf{x}, t)\]

In the pull scheme above, this translates to setting the incoming stream \(f_i\) equal to the current cell's pre-collision distribution along the opposite trajectory (\(f_{\text{opposite}(i)}\)) if the source cell is a solid node.

Maintaining Stability in GLSL

Relaxation Parameter Limits: The numerical scheme requires \(\tau > 0.5\). In single-precision floating-point (standard fp32 on GPUs), running \(\tau\) values near \(0.5\) (low viscosity) quickly triggers numerical instabilities and non-physical negative distributions. Typically, \(\tau \in [0.51, 1.2]\) preserves stability.
Mach Number Constraint: LBM operates under the weakly compressible Navier-Stokes assumption. Grid velocities must satisfy \(|\mathbf{u}| < 0.1 c_s\) to ensure negligible compressibility error and prevent divergence in the equilibrium distribution calculation.