Task: hotspot_2d
Category: Scientific simulation
Difficulty: 2
Tags: stencil,thermal,finite_difference,double_buffering,structured_grid
Source: Rodinia 3.1: https://github.com/HPC-FAIR/rodinia_3.1
Description: Two-dimensional thermal simulation from temperature and power grids using stencil updates.

Prompt template (YAML):
task_description: |
  Simulate transient temperature evolution on a 2D chip surface using the Rodinia Hotspot finite-difference model.

  Inputs are the initial temperature field `temp0[rows][cols]` and the power dissipation field `power[rows][cols]`.
  Run `iters` time steps and output the final temperature field.

interface: |
  ```c
  // Single end-to-end call. Do H2D + kernel + D2H and synchronize before returning.
  // Called multiple times (warmup + timed); must be idempotent.
  extern "C" void solution_compute(
      int rows,
      int cols,
      int iters,
      const float* temp0,
      const float* power,
      float* out_temp
  );
  ```

input_size_notes: |
  - `temp0` and `power` are row-major float arrays of length `rows*cols`.
  - **small**: 256x256, 50 iterations
  - **medium**: 768x768, 100 iterations
  - **large**: 1536x1536, 150 iterations

output_contract: |
  - Output the final temperature of every cell after `iters` iterations.
  - Output format is row-major float array written as one number per line.
  - Numerical tolerance: `atol=1e-3`, `rtol=1e-4`.

algorithm_background: |
  This is a structured-grid stencil update adapted from Rodinia Hotspot.
  Each iteration updates each cell using its current temperature, its immediate neighbors,
  the power input, and ambient-temperature dissipation.

  Interior cells use a 5-point stencil; edges and corners use one-sided boundary conditions.
  The computation naturally uses two buffers (`src` and `dst`) swapped after each iteration.

hints_l2: |
  - Keep `temp0` and `power` on GPU persistently after initialization.
  - Use a 2D thread block layout that maps naturally to the grid.
  - Tile the temperature grid into shared memory with halo cells.
  - Consider a separate fast path for interior tiles and a simpler path for boundaries.

hints_l1: |
  ### 1. Persistent device allocations
  Allocate device buffers for `src`, `dst`, and `power` inside `solution_compute`.

  ### 2. Iterative stencil kernel
  For each iteration, launch a 2D kernel that updates one cell per thread.

  ### 3. Boundary handling
  Rodinia Hotspot uses explicit formulas for corners and edges. Keep those conditions correct.

  ### 4. Shared-memory tiling
  Interior blocks can load a tile of `src` plus a one-cell halo to reduce global-memory traffic.

  ### 5. Final copy-back
  After the last iteration, copy the final buffer to `out_temp` on host.


CPU reference implementation (C):
// Adapted from the Rodinia Hotspot benchmark core solver
// Source repo: https://github.com/yuhc/gpu-rodinia
// Original file: openmp/hotspot/hotspot_openmp.cpp

#include <stdlib.h>
#include <string.h>

#ifdef __cplusplus
extern "C" {
#endif

typedef float FLOAT;

#define MAX_PD (3.0e6f)
#define PRECISION 0.001f
#define SPEC_HEAT_SI 1.75e6f
#define K_SI 100.0f
#define FACTOR_CHIP 0.5f

static const FLOAT t_chip = 0.0005f;
static const FLOAT chip_height = 0.016f;
static const FLOAT chip_width = 0.016f;
static const FLOAT amb_temp = 80.0f;

static int g_rows = 0;
static int g_cols = 0;
static int g_iters = 0;
static const FLOAT* g_temp0 = NULL;
static const FLOAT* g_power = NULL;
static FLOAT* g_buf_a = NULL;
static FLOAT* g_buf_b = NULL;

static void _orbench_old_init(int rows, int cols, int iters,
                   const float* temp0,
                   const float* power) {
    g_rows = rows;
    g_cols = cols;
    g_iters = iters;
    g_temp0 = temp0;
    g_power = power;

    size_t n = (size_t)rows * (size_t)cols;
    if (g_buf_a) free(g_buf_a);
    if (g_buf_b) free(g_buf_b);
    g_buf_a = (FLOAT*)malloc(n * sizeof(FLOAT));
    g_buf_b = (FLOAT*)malloc(n * sizeof(FLOAT));
}

static void _orbench_old_compute(float* out_temp) {
    if (!g_temp0 || !g_power || !g_buf_a || !g_buf_b || !out_temp ||
        g_rows <= 0 || g_cols <= 0 || g_iters <= 0) {
        return;
    }

    int row = g_rows;
    int col = g_cols;
    size_t n = (size_t)row * (size_t)col;
    memcpy(g_buf_a, g_temp0, n * sizeof(FLOAT));

    FLOAT grid_height = chip_height / (FLOAT)row;
    FLOAT grid_width  = chip_width  / (FLOAT)col;
    FLOAT Cap = FACTOR_CHIP * SPEC_HEAT_SI * t_chip * grid_width * grid_height;
    FLOAT Rx  = grid_width  / (2.0f * K_SI * t_chip * grid_height);
    FLOAT Ry  = grid_height / (2.0f * K_SI * t_chip * grid_width);
    FLOAT Rz  = t_chip / (K_SI * grid_height * grid_width);
    FLOAT max_slope = MAX_PD / (FACTOR_CHIP * t_chip * SPEC_HEAT_SI);
    FLOAT step = PRECISION / max_slope / 1000.0f;
    FLOAT Rx_1 = 1.0f / Rx;
    FLOAT Ry_1 = 1.0f / Ry;
    FLOAT Rz_1 = 1.0f / Rz;
    FLOAT Cap_1 = step / Cap;

    FLOAT* src = g_buf_a;
    FLOAT* dst = g_buf_b;

    for (int iter = 0; iter < g_iters; ++iter) {
        for (int r = 0; r < row; ++r) {
            for (int c = 0; c < col; ++c) {
                int idx = r * col + c;
                FLOAT delta;

                if (r == 0 && c == 0) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + 1] - src[idx]) * Rx_1
                        + (src[idx + col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (r == 0 && c == col - 1) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx - 1] - src[idx]) * Rx_1
                        + (src[idx + col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (r == row - 1 && c == col - 1) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx - 1] - src[idx]) * Rx_1
                        + (src[idx - col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (r == row - 1 && c == 0) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + 1] - src[idx]) * Rx_1
                        + (src[idx - col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (r == 0) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + 1] + src[idx - 1] - 2.0f * src[idx]) * Rx_1
                        + (src[idx + col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (r == row - 1) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + 1] + src[idx - 1] - 2.0f * src[idx]) * Rx_1
                        + (src[idx - col] - src[idx]) * Ry_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (c == 0) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + col] + src[idx - col] - 2.0f * src[idx]) * Ry_1
                        + (src[idx + 1] - src[idx]) * Rx_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else if (c == col - 1) {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + col] + src[idx - col] - 2.0f * src[idx]) * Ry_1
                        + (src[idx - 1] - src[idx]) * Rx_1
                        + (amb_temp - src[idx]) * Rz_1);
                } else {
                    delta = Cap_1 * (g_power[idx]
                        + (src[idx + col] + src[idx - col] - 2.0f * src[idx]) * Ry_1
                        + (src[idx + 1] + src[idx - 1] - 2.0f * src[idx]) * Rx_1
                        + (amb_temp - src[idx]) * Rz_1);
                }

                dst[idx] = src[idx] + delta;
            }
        }
        FLOAT* tmp = src;
        src = dst;
        dst = tmp;
    }

    memcpy(out_temp, src, n * sizeof(FLOAT));
}



#ifdef __cplusplus
}
#endif

// ── Unified compute_only wrapper (auto-migrated) ──
void solution_compute(int rows, int cols, int iters, const float* temp0, const float* power, float* out_temp) {
    _orbench_old_init(rows, cols, iters, temp0, power);
    _orbench_old_compute(out_temp);
}