Task: sph_position
Category: Scientific simulation
Difficulty: 1
Tags: sph,particle_update,embarrassingly_parallel,position_integration
Source: DualSPHysics: https://github.com/DualSPHysics/DualSPHysics
Description: Updates SPH particle positions and velocities for one integration step.

Prompt template (YAML):
# SPH Particle Position Update -- Prompt Template

task_description: |
  Update particle positions in a Smoothed Particle Hydrodynamics (SPH) simulation.
  Each particle has a 3D position (x, y, z) and a displacement vector (movx, movy, movz).
  For each particle: compute new_pos = old_pos + displacement, then compute the cell
  index from the new position using floor(pos / cell_size).

  This is the position integration step from DualSPHysics (Crespo et al. 2015),
  specifically the KerUpdatePos / KerComputeStepPos kernels. Each particle is
  completely independent, making this embarrassingly parallel.

  Based on: DualSPHysics JSphGpu_ker.cu, KerUpdatePos.
  Reference: Crespo et al. "DualSPHysics: Open-source parallel CFD solver based on
  SPH", Computer Physics Communications, 2015.

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 N,
      const float* posxy_x,
      const float* posxy_y,
      const float* posz,
      const float* movxy_x,
      const float* movxy_y,
      const float* movz,
      float cell_size,
      double* out_x,
      double* out_y,
      double* out_z,
      int* out_cell
  );
  ```

input_size_notes: |
  | Size   | N (particles) |
  |--------|--------------|
  | small  | 50,000       |
  | medium | 500,000      |
  | large  | 5,000,000    |

  - Each particle is independent -- one thread per particle
  - Positions are float64 (double precision, matching DualSPHysics)
  - Displacements are small (dt * velocity, typical SPH timestep ~1e-4)
  - Cell index is linearized as: floor(x/cell_size) + floor(y/cell_size)*1000 + floor(z/cell_size)*1000000

output_contract: |
  - out_x[i] = posxy_x[i] + movxy_x[i]
  - out_y[i] = posxy_y[i] + movxy_y[i]
  - out_z[i] = posz[i] + movz[i]
  - out_cell[i] = floor(out_x[i]/cell_size) + floor(out_y[i]/cell_size)*1000 + floor(out_z[i]/cell_size)*1000000
  - Tolerance: atol=1e-6, rtol=1e-5
  - solution_compute must be idempotent and synchronize before returning

algorithm_background: |
  **SPH Position Update (Symplectic Euler / Verlet)**

  In SPH fluid simulation, the position update step is:
  ```
  x_new = x_old + dt * v
  ```
  where dt*v is provided as the displacement vector (movx, movy, movz).

  After updating position, the cell index must be recomputed for the
  spatial data structure used in neighbor search:
  ```
  cx = floor(x / cell_size)
  cy = floor(y / cell_size)
  cz = floor(z / cell_size)
  cell = cx + cy * 1000 + cz * 1000000
  ```

  This kernel is memory-bandwidth bound on GPU: each particle reads 6 doubles
  (48 bytes) and writes 3 doubles + 1 int (28 bytes).

hints_l2: |
  - Each particle is independent -- one thread per particle, trivially parallel
  - All arrays are contiguous -- coalesced memory access is natural
  - Pre-upload all arrays to GPU inside solution_compute, copy results back in solution_compute
  - This is bandwidth-bound: minimize memory transactions
  - Consider using vectorized loads (double2) to match DualSPHysics data layout

hints_l1: |
  ### 1. Data Transfer Strategy
  In solution_compute setup:
  - Allocate device buffers for 6 input arrays (double) + 3 output arrays (double) + 1 output (int)
  - cudaMemcpy all inputs to device (done once)

  In solution_compute:
  - Launch kernel, then cudaMemcpy results back to host

  ### 2. Kernel Structure
  ```cuda
  __global__ void KerUpdatePos(int N, const double* px, const double* py,
      const double* pz, const double* mx, const double* my, const double* mz,
      double cell_size, double* ox, double* oy, double* oz, int* ocell) {
      int i = blockIdx.x * blockDim.x + threadIdx.x;
      if (i >= N) return;
      double rx = px[i] + mx[i];
      double ry = py[i] + my[i];
      double rz = pz[i] + mz[i];
      ox[i] = rx; oy[i] = ry; oz[i] = rz;
      ocell[i] = (int)floor(rx/cell_size) + (int)floor(ry/cell_size)*1000
               + (int)floor(rz/cell_size)*1000000;
  }
  ```

  ### 3. Launch Configuration
  - 256 threads per block
  - Grid: (N + 255) / 256 blocks
  - No shared memory needed


CPU reference implementation (C):
// cpu_reference.c -- SPH particle position update (CPU baseline)
//
// Faithfully ported from DualSPHysics JSphGpu_ker.cu
//   KerUpdatePos / KerComputeStepPos
//
// Simplified: no periodic boundaries, no floating bodies.
// For each particle: new_pos = old_pos + displacement, recompute cell index.
// NO file I/O, NO main(). All I/O handled by task_io_cpu.c.
//
// Reference: Crespo et al. "DualSPHysics: Open-source parallel CFD solver
// based on SPH", Computer Physics Communications, 2015.

#include <math.h>
#include <stdlib.h>

// ===== Module-level static state =====
static int    g_N;
static const float* g_pos_x;
static const float* g_pos_y;
static const float* g_pos_z;
static const float* g_mov_x;
static const float* g_mov_y;
static const float* g_mov_z;
static float g_cell_size;

// ===== Ported from DualSPHysics KerUpdatePos (simplified, non-periodic) =====
// Original: template<bool periactive> __device__ void KerUpdatePos(...)
// Applies displacement to particle position and computes new cell index.

static void KerUpdatePos(double px, double py, double pz,
                         double movx, double movy, double movz,
                         float cell_size,
                         double* out_x, double* out_y, double* out_z,
                         int* out_cell_x, int* out_cell_y, int* out_cell_z)
{
    /* Apply displacement (matches: rpos.x+=movx; rpos.y+=movy; rpos.z+=movz;) */
    double rx = px + movx;
    double ry = py + movy;
    double rz = pz + movz;

    *out_x = rx;
    *out_y = ry;
    *out_z = rz;

    /* Compute cell indices (matches dcell computation in DualSPHysics) */
    *out_cell_x = (int)floor(rx / cell_size);
    *out_cell_y = (int)floor(ry / cell_size);
    *out_cell_z = (int)floor(rz / cell_size);
}

// ===== Public interface =====

static void _orbench_old_init(int N,
                   const float* posxy_x, const float* posxy_y,
                   const float* posz,
                   const float* movxy_x, const float* movxy_y,
                   const float* movz,
                   float cell_size)
{
    g_N = N;
    g_pos_x = posxy_x;
    g_pos_y = posxy_y;
    g_pos_z = posz;
    g_mov_x = movxy_x;
    g_mov_y = movxy_y;
    g_mov_z = movz;
    g_cell_size = cell_size;
}

static void _orbench_old_compute(int N,
                      double* out_x, double* out_y, double* out_z,
                      int* out_cell)
{
    int i;
    for (i = 0; i < N; i++) {
        int cx, cy, cz;
        KerUpdatePos(g_pos_x[i], g_pos_y[i], g_pos_z[i],
                     g_mov_x[i], g_mov_y[i], g_mov_z[i],
                     g_cell_size,
                     &out_x[i], &out_y[i], &out_z[i],
                     &cx, &cy, &cz);
        /* Encode cell as linear index: cx + cy*1000 + cz*1000000 */
        out_cell[i] = cx + cy * 1000 + cz * 1000000;
    }
}

// ── Unified compute_only wrapper (auto-migrated) ──
void solution_compute(int N, const float* posxy_x, const float* posxy_y, const float* posz, const float* movxy_x, const float* movxy_y, const float* movz, float cell_size, double* out_x, double* out_y, double* out_z, int* out_cell) {
    _orbench_old_init(N, posxy_x, posxy_y, posz, movxy_x, movxy_y, movz, cell_size);
    _orbench_old_compute(N, out_x, out_y, out_z, out_cell);
}