GLSL Smoke Buoyancy: Coupling Temperature and Density

This article explains how to simulate realistic smoke buoyancy in real-time fluid solvers by coupling temperature and density fields inside a GLSL compute shader. It covers the governing physical approximation, the mathematical formulation for buoyant acceleration, and the step-by-step implementation of an external force pass within a modern GPU compute pipeline.

The Physics: The Boussinesq Approximation

In real-time smoke simulation based on the incompressible Navier-Stokes equations, buoyancy is typically modeled using the Boussinesq approximation. This approximation assumes that variations in fluid density are negligible except where they contribute to the gravitational buoyant force.

Warm smoke rises because thermal expansion makes it lighter than surrounding air, while heavy particulate concentrations (smoke density) pull it downward under gravity. The upward buoyant force \(F_{\text{buoy}}\) acting along the vertical axis is defined as:

\[F_{\text{buoy}} = \left( -\alpha \rho + \beta (T - T_{\text{amb}}) \right) \hat{y}\]

  • \(\rho\): The local smoke density.
  • \(T\): The local temperature.
  • \(T_{\text{amb}}\): The ambient room temperature.
  • \(\alpha\): A positive scalar weighting the downward gravitational pull of dense smoke.
  • \(\beta\): A positive scalar controlling thermal lift (thermal expansion coefficient).
  • \(\hat{y}\): The normalized vertical direction vector (typically \((0, 1, 0)\)).

When the gas is hot (\(T > T_{\text{amb}}\)), the thermal term dominates and accelerates the fluid upward. As the smoke cools or accumulates heavy density, the downward pull balances or overcomes the thermal lift, creating natural swirling plumes and mushroom caps.

Storing Fluid States on the GPU

In a typical GLSL pipeline, the 3D fluid domain is represented using 3D textures. You will generally use:

  • velocityTexture (format RGBA16F or RGBA32F): Stores the 3D velocity vector \((u, v, w)\) in the RGB channels.
  • densityTexture (format R16F or R32F): Stores scalar smoke concentration \(\rho\).
  • temperatureTexture (format R16F or R32F): Stores scalar temperature \(T\).

Alternatively, density and temperature can be packed into a single two-channel RG16F texture to reduce memory bandwidth and cache misses during texture sampling.

GLSL Compute Shader Implementation

The buoyancy calculation is executed as an explicit force accumulation step. It runs after the advection stage and before the pressure projection (divergence-free projection) stage of the fluid solver.

#version 450 core

layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;

// Textures
layout(binding = 0, rgba16f) uniform readonly image3D u_VelocityIn;
layout(binding = 1, r16f)    uniform readonly image3D u_Density;
layout(binding = 2, r16f)    uniform readonly image3D u_Temperature;
layout(binding = 3, rgba16f) uniform writeonly image3D u_VelocityOut;

// Simulation parameters
uniform float u_DeltaTime;
uniform float u_AmbientTemperature;
uniform float u_Alpha; // Density weight
uniform float u_Beta;  // Temperature buoyancy weight

void main()
{
    ivec3 cell = ivec3(gl_GlobalInvocationID.xyz);
    ivec3 size = imageSize(u_VelocityIn);

    // Boundary check
    if (cell.x >= size.x || cell.y >= size.y || cell.z >= size.z) {
        return;
    }

    // Read current state
    vec3 velocity = imageLoad(u_VelocityIn, cell).xyz;
    float density = imageLoad(u_Density, cell).r;
    float temp    = imageLoad(u_Temperature, cell).r;

    // Calculate vertical buoyant force
    float thermalLift = u_Beta * (temp - u_AmbientTemperature);
    float densityWeight = u_Alpha * density;
    float verticalForce = thermalLift - densityWeight;

    // Apply explicit Euler integration along the vertical axis
    velocity.y += verticalForce * u_DeltaTime;

    // Output updated velocity field
    imageStore(u_VelocityOut, cell, vec4(velocity, 0.0));
}

Integrating into the Simulation Loop

To maintain stability and realistic motion, integrate the buoyancy compute pass into the larger simulation pipeline in the following sequence:

  1. Advection: Advect velocity, density, and temperature using semi-Lagrangian backtracing or higher-order schemes (e.g., MacCormack).
  2. Thermal Cooling and Diffusion: Apply heat dissipation to gradually return temperatures toward \(T_{\text{amb}}\), along with optional thermal and density diffusion.
  3. Buoyancy Force Application: Dispatch the compute shader shown above to update the velocity field based on current temperature and density values.
  4. Pressure Projection (Poisson Solve): Compute the divergence of the updated velocity field and resolve the pressure Poisson equation (using Jacobi iterations, Gauss-Seidel, or Multigrid). Subtract the pressure gradient from the velocity to enforce mass conservation (incompressibility).

Placing buoyancy right before the pressure solve is critical. The pressure projection step removes the divergence introduced by the localized upward acceleration, converting pure vertical velocity into the circular vorticity patterns characteristic of rising smoke plumes.