# The OpenCL (GPU) implementation

This page explains how the PIV algorithm is parallelised on the GPU in
`dpivsoft/Cl_DPIV.py`. The algorithm itself — normalized FFT cross-correlation,
sub-pixel peak fitting, multi-pass window deformation — is identical to the CPU
implementation and is described on the *PIV algorithm* page; here the subject
is the execution model, the kernels, and the memory management.

The implementation is **vendor-neutral**: OpenCL through `pyopencl`, with
[reikna](http://reikna.publicfields.net/) providing the batched GPU FFT and the
device abstraction. It runs on any GPU (or CPU device) with an OpenCL runtime —
no CUDA dependency.

---

## 1. Why PIV maps well to a GPU

The expensive part of PIV is embarrassingly parallel at two levels: **across
windows** (thousands of interrogation windows, each correlated independently)
and **across pixels** (slicing, mean subtraction, normalization, deformation
and the FFT butterflies are per-pixel operations).

The CPU implementation loops over windows one at a time; the GPU implementation
instead stores **all windows of a pass as one 3-D array** of shape

$$
(N_{\text{boxes}},\; B_y,\; B_x)
$$

and every pipeline stage processes the whole stack in a single kernel launch —
the batched FFT transforms axes $(1, 2)$, all windows at once. This
restructuring, not any algorithmic change, is where the speed-up (order
$\times 100$ over the pure-Python path) comes from.

---

## 2. Software stack

| layer | role |
|---|---|
| `pyopencl` | OpenCL context, buffers, kernel launch |
| `reikna` | device thread abstraction, batched FFT / FFTShift, temp-array manager |
| `GPU_Kernels/*.cl` | hand-written OpenCL C kernels (one file per stage) |
| `Classes.GPU` | namespace holding all compiled kernels and device arrays |

All arithmetic is single precision (`float32` / `complex64`) — the native fast
path on GPUs. Results agree with the double-precision CPU implementation within
tolerances verified by `tests/test_GPUconsistency.py`, which compares grid
generation, slicing, normalization, correlation, peak finding, derivatives,
interpolation and image deformation between the two implementations.

---

## 3. The kernel inventory

Each `.cl` file implements one pipeline stage; its CPU counterpart is the
function of the same role in `DPIV.py`:

| kernel | stage | CPU counterpart |
|---|---|---|
| `Slice.cl` | cut all sub-windows into the 3-D stack | array slicing |
| `SubMean.cl` | per-window mean (tree reduction, local memory) | `np.sum` per window |
| `Normalize_Img.cl` | subtract mean, divide by intensity norm | mean-subtraction + $\sigma$ |
| reikna `FFT` / `FFTShift` | batched 2-D FFT over all windows | `np.fft.fft2` / `fftshift` |
| `conjugate_inplace.cl` | conjugate the first spectrum in place | `np.conj` |
| `multiply_them.cl` | spectrum product (also reused, with `float32×bool`, as the masking kernel) | spectrum product |
| `find_peak.cl` | two-peak search + Westerweel correction + Gaussian sub-pixel fit | `find_peaks` |
| `median_filter.cl` | median-test outlier replacement | `median_filter` |
| `Jacobian.cl` + `box_blur.cl` | velocity gradients + 3×3 smoothing | `jacobian_matrix` |
| `deform_image.cl` | symmetric window deformation (bilinear) | `deform_image` |
| `interpolation.cl` | pass-1 field → pass-2 grid | `interpolations` |
| `gaussian_filter.cl` | first-pass image smoothing | `gaussian_filter` |
| `Weighting.cl` | window weighting | `weight_function` |
| `directCorrelation.cl` + `find_peak_direct.cl` | direct-correlation first pass | `corrDirect1` |
| `ini_index.cl` | zero the per-pixel displacement buffers | array initialisation |

---

## 4. Set-up: three calls, once each

```python
import dpivsoft.Cl_DPIV as Cl_DPIV
from dpivsoft.Classes import Parameters, GPU

thr = Cl_DPIV.select_Platform(0)        # or "selection" for an interactive list
Cl_DPIV.compile_Kernels(thr)            # compile geometry-independent kernels
Cl_DPIV.initialization(width, height, thr)   # allocate buffers, compile the rest
```

- **{func}`~dpivsoft.Cl_DPIV.select_Platform`** creates the reikna `Thread`
  on the chosen OpenCL platform/device.
- **{func}`~dpivsoft.Cl_DPIV.compile_Kernels`** compiles every kernel whose
  source does not depend on the PIV geometry. Once per device.
- **{func}`~dpivsoft.Cl_DPIV.initialization`** does everything geometry-dependent: builds and uploads
  the PIV grid and parameter blocks, allocates all device arrays, compiles the
  two remaining kernels ([§5](#ocl-sec-5)), and plans the batched FFTs. Call it
  again whenever the image size or the PIV parameters change.

(ocl-sec-5)=
## 5. Geometry-dependent kernels and work-group limits

Two kernels are compiled inside `initialization` rather than `compile_Kernels`,
because their **local (shared) memory arrays need compile-time sizes** that
depend on the window geometry (the sources are Mako templates):

- **`SubMean`** computes each window's mean with a binary-tree reduction in
  local memory: the work-group size must be a power of two dividing the window
  pixel count.
- **`find_peak`** searches the correlation maximum with the same reduction
  pattern over the search window, then applies the Westerweel correction and
  the three-point Gaussian fit in-kernel.

Both are compiled through a clamping loop: the requested work-group size is the
largest power of two ≤ min(window pixels, device maximum), but drivers may
report a **smaller per-kernel limit** (`CL_KERNEL_WORK_GROUP_SIZE`, bound by
register and local-memory usage). The loop compiles, reads the kernel's real
limit, and recompiles clamped to it if needed — a no-op on lenient drivers; on
stricter runtimes (e.g. Mesa's rusticl) it is what makes the launch valid.

---

## 6. Memory management

Three decisions keep device memory small and allocation-free at run time:

- **Everything is allocated once**, in `initialization`; the per-frame
  `processing` call performs zero allocations.
- **Temporary arrays share memory** via reikna's `ZeroOffsetManager` with
  explicit dependency lists: two temp arrays may occupy the same physical
  memory unless one depends on the other being alive, so pass-1 and pass-2
  intermediates overlap. One visible consequence: buffers can hold **stale
  data** from the other pass, so the per-pixel displacement buffers are
  explicitly zeroed (`ini_index`) at the start of each frame.
- **The image buffers are persistent**: `GPU.img1` / `GPU.img2` are allocated
  once and overwritten with `.set()` for every new pair.

## 7. Hiding the file I/O

`processing(img1_name, img2_name, thr)` takes the file names of the **next**
image pair, not the current one. Kernel launches are asynchronous: after the
last kernel of the current pair is enqueued, the host thread — while the GPU
queue drains — loads the next pair from disk, then synchronises and uploads it.
Disk I/O is thereby overlapped with GPU compute, which matters for long image
sequences:

```python
GPU.img1.set(Img1)            # upload the FIRST pair manually
GPU.img2.set(Img2)

for i in range(0, len(files), 2):
    nxt1, nxt2 = next_pair_names(i)          # names of the NEXT pair
    Cl_DPIV.processing(nxt1, nxt2, thr)      # compute current + preload next
    x = GPU.x2.get(); y = GPU.y2.get()
    u = GPU.u2_f.get(); v = GPU.v2_f.get()   # final (median-filtered) field
```

Results stay in device memory as attributes of the `GPU` class; `.get()`
transfers them to NumPy arrays.

---

## 8. The pipeline, stage by stage

`processing` mirrors the CPU control flow exactly. Per frame:

1. Zero the displacement-index buffers (`ini_index`); mask and/or Gaussian-blur
   the images if enabled.
2. **Pass 1** (repeated `no_iter_1` times): slice sub-windows → (from the
   second iteration: median filter → Jacobian + box blur → deform) → per-window
   mean + normalize → batched FFT → conjugate → spectrum product → inverse FFT
   → FFT-shift → `find_peak`. With `direct_calc`, the first iteration instead
   runs `directCorrelation` + `find_peak_direct`.
3. Interpolate the pass-1 field onto the fine grid (`interpolation`).
4. **Pass 2** (repeated `no_iter_2` times): Jacobian + box blur → deform →
   normalize → FFT correlation → `find_peak` → median filter (+ mask).
5. Load the next image pair from disk, synchronise, upload.

One deliberate difference from the CPU path: window **weighting is disabled**
on the GPU (a known bug in the weighting kernel); if requested, `processing`
turns it off and prints a warning.

---

## 9. Function reference

| function | role |
|----------|------|
| {func}`~dpivsoft.Cl_DPIV.select_Platform` | choose the OpenCL platform/device, return the reikna thread |
| {func}`~dpivsoft.Cl_DPIV.compile_Kernels` | compile all geometry-independent kernels (once per device) |
| {func}`~dpivsoft.Cl_DPIV.initialization` | grid + buffers + geometry-dependent kernels + FFT plans (once per configuration) |
| {func}`~dpivsoft.Cl_DPIV.processing` | run the full two-pass pipeline on the pair in GPU memory; preload the named next pair |

---

## References

1. J. Aguilar-Cabello, L. Parras, C. del Pino (2022), "[DPIVSoft-OpenCL: A
   multicore CPU–GPU accelerated open-source code for 2D Particle Image
   Velocimetry](https://www.sciencedirect.com/science/article/pii/S2352711022001741),"
   *SoftwareX* **20**, 101256,
   [doi:10.1016/j.softx.2022.101256](https://doi.org/10.1016/j.softx.2022.101256).
   — The code paper: this OpenCL parallelisation, its design and its
   validation.
2. P. Meunier, T. Leweke (2003), "[Analysis and treatment of errors due to high
   velocity gradients in particle image
   velocimetry](https://doi.org/10.1007/s00348-003-0673-2)," *Exp. Fluids*
   **35**, 408–421. — The underlying algorithm (see the *PIV algorithm* page).
3. [Khronos OpenCL specification](https://registry.khronos.org/OpenCL/) — the
   portable compute standard this implementation targets through `pyopencl`.