Building a Real-Time Fluid Simulator on the GPU (CUDA + OpenGL) — Heat Simulation (Part 4)

If you have not read Part 3 of the series, please read the article first, as it presents the project setup upon which we will build in this chapter.

Please refer to the implementation on GitHub for the full code for this chapter.

Before we jump into Navier–Stokes and full fluid dynamics, it’s useful to warm up (pun intended) with a simpler physical system: heat diffusion.

The heat equation is one of the friendliest PDEs to simulate on a grid, and it gives us the perfect opportunity to introduce concepts that will later become essential: the Laplacian, numerical time stepping, buffer swapping, stability constraints, and visualizing scalar fields in CUDA.

Let’s get started.

1. The Physics: The Heat Equation

The heat diffusion equation in 2D is:

Where:

• T(x, y, t) is the temperature field
• α is the diffusion coefficient
• ∇²T is the Laplacian (a measure of how much the local value differs from its neighbors)

Intuition:

• If a cell is warmer than its neighbors → it cools down.
• If it’s colder → it warms up.
• Everything slowly smooths out.

This is exactly the kind of system that GPUs are great at: local, grid-based, massively parallel updates.

2. Discretizing the Laplacian

On a grid, we approximate the Laplacian with the 5-point stencil:

In CUDA, we compute this using simple neighbor accesses. This is a local operation → one CUDA thread computes one cell.

3. The CUDA Kernel: Computing One Time Step

Here’s the core of the heat simulation:

// Compute the discrete Laplacian for a 2D grid at (x, y)
// The Laplacian approximates the second spatial derivative, ∇²u, which governs diffusion.
// Each neighbor contributes the difference from the center cell.
__device__ __forceinline__ float laplacian(const float* grid, int x, int y, int width, int height) {
    float center = grid[y * width + x];
    float sum = 0.0f;

    if (x > 0) sum += grid[y * width + (x - 1)] - center;
    if (x < width - 1) sum += grid[y * width + (x + 1)] - center;
    if (y > 0) sum += grid[(y - 1) * width + x] - center;
    if (y < height - 1) sum += grid[(y + 1) * width + x] - center;

    return sum;
}

// Advance one timestep of the heat simulation using explicit Euler method:
// newValue = oldValue + dt * diffusion * Laplacian(oldValue)
// This models heat spreading from each cell to neighbors.
__global__ void heatStepKernel(const float* current, float* next, int width, int height, float dt, float diffusion, float sourceValue) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;

    int idx = y * width + x;
    float lap = laplacian(current, x, y, width, height);

    // Update value according to diffusion PDE: u_t = D * ∇²u
    next[idx] = current[idx] + diffusion * lap * dt;

    // Apply constant heat source at left boundary (x == 0)
    // This ensures the simulation always has a “hot edge”
    if (x == 0) {
        next[idx] = sourceValue;
    }
}

All we did was translate the Laplacian approximation to code in the laplacian kernel, and use it the Euler update step in the heatStepKernel . Note that I added boundary checks to make sure we do not overflow. Also note that I am simulating a heat source on the left by constantly pumping new heat in the leftmost cells.

This is why we get the following pattern:

Heat diffusing towards the right from a heat source on the left.

Note that because each cell needs information from its neighboring cells, if we were to update in-place, the evolution of the simulation would depend on which threads finished their work first and yield wrong results. Instead, we have to have to write to a second buffer to avoid this pitfall. In practice, we maintain two grids:

• m_devCurrent — the state at time t
• m_devNext — the state at time t + Δt: this is where we write the update every step

After each step, we swap the buffers so that we do not have to copy the data or have complex bookkeeping logic:

std::swap(m_devCurrent, m_devNext);

This is the classic double-buffering pattern used in many simulations: one buffer stores the “current” state while the next state is written into another buffer. After each step, the two are swapped. Your GPU uses a very similar idea when rendering frames — drawing into a back buffer while the fully-rendered front buffer is being displayed — so the monitor never shows a half-updated image.

4. Mapping Heat Values to Colors

Heat itself is just a scalar field, so we convert it to an uchar4 RGBA pixel:

// Clamp values to lie between 0 and 1
__device__ __forceinline__ float saturatef(float x) {
    return fminf(fmaxf(x, 0.0f), 1.0f);
}

// Simple linear mapping: low values -> blue, high values -> red
__device__ __forceinline__ float4 heatToColor(float value) {
    // Saturate value to [0,1] to avoid color overflow
    return make_float4(saturatef(value),  saturatef(value*0.3f), saturatef(1.0f - value), 1.0f);
}

This gives us a simple blue → red colormap.

Then a second kernel writes the pixel data directly into the CUDA-OpenGL shared PBO:

__global__ void heatToColorKernel(
    const float* heat,
    uchar4* buffer,
    int width,
    int height)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;

    int idx = y * width + x;
    buffer[idx] = floatToUchar4(heatToColor(heat[idx]));
}

5. The Full Step Function (Putting It Together)

Finally, we can have a look at what the whole simulation loop looks like. Two kernels, one after the other:

1. Update the temperature field
2. Convert temperature to pixels in the PBO

OpenGL will then upload that PBO to a texture and draw it.

void HeatSimulation::step(uchar4* pbo) {
    // CUDA kernel launch dimensions: 16x16 threads per block
    dim3 block(16, 16);
    dim3 grid((m_width + block.x - 1) / block.x, (m_height + block.y - 1) / block.y);

    // 1) Advance simulation one timestep
    heatStepKernel<<<grid, block>>>(m_devCurrent, m_devNext, m_width, m_height, m_dt, m_diffusion, m_sourceHeat);
    cudaDeviceSynchronize();

    // Check for errors
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        std::cerr << "CUDA post-sync error: " << cudaGetErrorString(err) << std::endl;
    }

    // 2) Convert heat values to RGBA colors for rendering
    heatToColorKernel<<<grid, block>>>(m_devNext, pbo, m_width, m_height);
    cudaDeviceSynchronize();

    err = cudaGetLastError();
    if (err != cudaSuccess) {
        std::cerr << "CUDA post-sync error: " << cudaGetErrorString(err) << std::endl;
    }

    // Swap grids for next timestep (ping-pong buffer)
    std::swap(m_devCurrent, m_devNext);
}

Note that it’s good practice to always check for the last potential CUDA error as it helps with debugging when things go wrong. Running the heat simulation should now show:

• A red band on the left injecting heat
• Warmth diffusing smoothly across the grid
• A nice gradient that evolves over time

6. A Note on Numerical Stability

The explicit Euler update is simple, but it comes with a numerical restriction:

Meaning:

• Too large a timestep → unstable simulation
• Too high diffusion → the temperature explodes

In my case, I chose:

• dt = 0.6f
• diffusion = 0.2f

The grid spacing dx can be considered equal to 1, plugging in the values, we get:

which clearly holds.

In a future article, we will explore more stable methods such as “Implicit Euler”, “Semi-Implicit methods”, “ADI schmes” etc. but these are more complex to implement.

Coming Up Next

In the next article, we will finally be setting up the long-due fluid simulation. Similarly to this section, there will be quite a bit of math and even more kernels! How exciting — hopefully.

Let’s keep going!