D3Q19 LBM Stream and Collide in GLSL

This article provides an overview of implementing the stream-and-collide algorithm for the D3Q19 Lattice Boltzmann Method (LBM) in OpenGL Shading Language (GLSL). It explains memory organization using 3D textures or buffers, the GPU-friendly "pull" streaming convention, the reduction of populations into macroscopic density and velocity, and the local relaxation toward equilibrium using the Bhatnagar-Gross-Krook (BGK) collision operator.

Memory Representation and Double Buffering

The D3Q19 lattice model tracks 19 discrete particle velocity distributions (\(f_0\) through \(f_{18}\)) per grid node. In GLSL, this data is typically mapped across five RGBA32F 3D textures (providing 20 channels, leaving one unused) or stored as a flat array within a Shader Storage Buffer Object (SSBO) inside a compute shader.

To prevent read-write race conditions across adjacent threads, a double-buffering (ping-pong) scheme is required. The current time step reads exclusively from the read texture/buffer \(t\), and writes newly computed values to the write texture/buffer \(t + \Delta t\).

The "Pull" Streaming Scheme

While theoretical LBM describes streaming as "pushing" distributions to adjacent nodes, GPU parallelization requires a "pull" scheme. Each thread corresponds to a target lattice cell at coordinate \(\vec{x}\) and fetches the incoming populations that were directed toward it from neighboring cells at \(\vec{x} - \vec{c}_i\), where \(\vec{c}_i\) is the discrete lattice velocity vector:

Rest velocity (1 vector): \(\vec{c}_0 = (0, 0, 0)\)
Face neighbors (6 vectors): \(\vec{c}_{1..6} = (\pm 1, 0, 0), (0, \pm 1, 0), (0, 0, \pm 1)\)
Edge neighbors (12 vectors): \(\vec{c}_{7..18} = (\pm 1, \pm 1, 0), (\pm 1, 0, \pm 1), (0, \pm 1, \pm 1)\)

By sampling the neighbor coordinate \(\vec{x} - \vec{c}_i\) from the input state, streaming is performed implicitly during the read phase.

Calculating Macroscopic Moments

After pulling the 19 distribution values \(f_i\), the local fluid density \(\rho\) and macroscopic velocity \(\vec{u}\) are derived via summation:

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

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

In GLSL, this is carried out using vector dot products and additions over the cached local variables.

The BGK Collision Step

The standard single-relaxation-time Bhatnagar-Gross-Krook (BGK) model relaxes the distribution functions toward an equilibrium state \(f_i^{eq}\). The equilibrium distribution is defined as:

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

For the D3Q19 lattice, the speed of sound constant is \(c_s^2 = \frac{1}{3}\), making \(c_s^4 = \frac{1}{9}\). The directional lattice weights \(w_i\) are:

\(w_0 = \frac{1}{3}\) (rest particle)
\(w_{1..6} = \frac{1}{18}\) (face connections, \(|\vec{c}_i|^2 = 1\))
\(w_{7..18} = \frac{1}{36}\) (edge connections, \(|\vec{c}_i|^2 = 2\))

The post-collision value \(f_i^*\) to be saved into the write buffer is calculated by relaxing \(f_i\) by a relaxation parameter \(\tau\):

\[f_i^* = f_i - \frac{1}{\tau} (f_i - f_i^{eq})\]

Here, \(\tau\) directly dictates the fluid's kinematic shear viscosity \(\nu = c_s^2 (\tau - 0.5) \Delta t\).

GLSL Implementation Routine

In a GLSL compute shader, the combined step is structured as follows:

#version 450
layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;

// 19 velocity vectors and weights defined as constants
const ivec3 c[19] = ivec3[19](...);
const float w[19] = float[19](...);

layout(binding = 0) uniform sampler3D u_PopulationsIn[5];
layout(binding = 1, rgba32f) uniform writeonly image3D u_PopulationsOut[5];

uniform float u_Tau;

void main() {
    ivec3 pos = ivec3(gl_GlobalInvocationID.xyz);

    float f[19];
    float rho = 0.0;
    vec3 momentum = vec3(0.0);

    // 1. Pull Streaming & Moment Accumulation
    for (int i = 0; i < 19; ++i) {
        ivec3 sourcePos = pos - c[i];
        
        // Unpack f[i] from the appropriate texture and channel
        f[i] = fetchPopulation(u_PopulationsIn, sourcePos, i);
        
        rho += f[i];
        momentum += f[i] * vec3(c[i]);
    }

    vec3 u = momentum / rho;
    float u_sq = dot(u, u);

    // 2. Collision & Write-back
    for (int i = 0; i < 19; ++i) {
        vec3 ci = vec3(c[i]);
        float ci_u = dot(ci, u);
        
        // Equilibrium computation (cs^2 = 1/3)
        float feq = w[i] * rho * (1.0 + 3.0 * ci_u + 4.5 * (ci_u * ci_u) - 1.5 * u_sq);
        
        // BGK Relaxation
        float f_star = f[i] - (f[i] - feq) / u_Tau;

        // Store result to output texture bindings
        storePopulation(u_PopulationsOut, pos, i, f_star);
    }
}

Handling boundary conditions, such as standard half-way bounce-back for no-slip solid obstacles, involves checking a voxel flag before collision. If a node is tagged as solid, the incoming streamed populations are reflected along their inverted direction vectors (\(\vec{c}_j = -\vec{c}_i\)) and written directly back into the output buffers without applying the BGK relaxation.