Simulating Radiative Heat Transfer with GLSL Compute Shaders

This article explores the end-to-end methodology for simulating radiative heat transfer across discretized surface meshes using OpenGL Shading Language (GLSL) compute shaders. It covers the mathematical formulation of radiative exchange based on enclosure theory, the translation of surface mesh data into GPU-friendly buffers, parallel view factor calculation using ray-tracing techniques, and the iterative numerical solvers used to compute steady-state or transient temperature fields entirely on the GPU.

Mathematical Formulation of Surface Radiation

Radiative exchange between opaque, diffuse-gray surfaces is governed by enclosure theory and the Stefan-Boltzmann law. The surface mesh is discretized into \(N\) planar elements (triangles or quadrilaterals). For each element \(i\), the radiosity \(J_i\) (total radiative energy leaving the surface per unit area) is defined as:

\[J_i = \epsilon_i \sigma T_i^4 + (1 - \epsilon_i) \sum_{j=1}^{N} F_{i \to j} J_j\]

where \(\epsilon_i\) is the surface emissivity, \(\sigma\) is the Stefan-Boltzmann constant, \(T_i\) is the absolute temperature, and \(F_{i \to j}\) is the geometric view factor describing the fraction of diffuse radiation leaving element \(i\) that directly intercepts element \(j\).

The net radiative heat flux \(q_i\) entering or leaving element \(i\) with surface area \(A_i\) is:

\[q_i = \frac{\epsilon_i}{1 - \epsilon_i} (\sigma T_i^4 - J_i) = \sum_{j=1}^N F_{i \to j}(J_i - J_j)\]

GPU Data Representation

Simulating this process in GLSL compute shaders requires organizing mesh geometry and thermodynamic properties into Shader Storage Buffer Objects (SSBOs).

  1. Element Geometry SSBO: Stores element centroids, surface normal vectors, and surface areas (\(A_i\)).
  2. Thermal State SSBOs: Two ping-pong buffers containing element temperatures (\(T_i\)), radiosities (\(J_i\)), and thermophysical properties (emissivity \(\epsilon_i\), heat capacity \(C_i\)).
  3. Spatial Acceleration SSBO: A flattened Bounding Volume Hierarchy (BVH) or uniform spatial grid used to resolve geometric occlusion during view factor computation.

View Factor Computation in GLSL

The most computationally demanding phase is calculating the view factor matrix \(F_{ij}\), which requires solving visibility and geometric configuration:

\[F_{i \to j} = \frac{1}{A_i} \int_{A_i} \int_{A_j} \frac{\cos\theta_i \cos\theta_j}{\pi r^2} V_{ij} \, dA_j \, dA_i\]

where \(r\) is the distance between elements, \(\theta_i\) and \(\theta_j\) are the angles between the surface normals and the connecting vector, and \(V_{ij} \in \{0, 1\}\) is a binary visibility function.

Compute shaders evaluate this using one of two primary GPU-accelerated approaches:

  • Ray-Traced Visibility: For moderate mesh resolutions, each compute shader thread corresponds to a pair \((i, j)\) or evaluates all \(j\) elements for a single element \(i\). Threads calculate the differential view factor based on centroid distances and evaluate \(V_{ij}\) by traversing the BVH stored in the spatial acceleration SSBO to detect intersecting occluders.
  • Hemicube Projection: For complex scenes where explicit pair-checking is too costly, each thread maps a virtual hemicube over an element's centroid. Surrounding primitives are projected onto the hemicube faces to aggregate incoming radiation weights without requiring explicit \(N \times N\) matrix storage.

Numerical Iteration via Compute Passes

Once view factors (or direct radiative exchanges) are resolved, the nonlinear radiosity equations must be solved iteratively. This is structured as a sequence of compute shader dispatches:

  1. Radiosity Evaluation Pass: A compute shader calculates new radiosity values \(J_i\) from current temperatures using a Jacobi iteration method. Because each element reads from the current read-state SSBO and writes to the write-state SSBO, race conditions are avoided across threads.
  2. Transient Heat Equation Pass: If simulating time-dependent heat transfer, the net heat flux \(q_i\) updates the element's temperature over a discrete time step \(\Delta t\):

\[T_i^{t+\Delta t} = T_i^t + \frac{q_i A_i}{C_i} \Delta t\]

  1. Ping-Pong Buffer Swap: The application swaps the bindings for the read and write SSBOs, repeating the dispatch loop until convergence criteria (residual radiative tolerance) or the target simulation time is met.

Performance Optimizations

To handle fine surface meshes efficiently, modern GLSL implementations utilize:

  • Subgroup Intrinsics and Shared Memory: Groups of threads process nearby surface patches cooperatively, loading patch properties into shared memory to minimize global memory bandwidth saturation.
  • Cut-Off Thresholds and Clustering: Distant surface elements are grouped into low-resolution thermal clusters, and element pairs where the distance \(r\) renders \(F_{ij} \approx 0\) are culled early before executing BVH intersection tests.
  • Hardware Ray Tracing: In modern pipelines supporting GLSL ray query extensions (GL_EXT_ray_query), compute shaders utilize dedicated RT-cores for occlusion tests (\(V_{ij}\)), accelerating visibility queries compared to software-based BVH traversals.