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 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

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

  • select_Platform() creates the reikna Thread on the chosen OpenCL platform/device.

  • compile_Kernels() compiles every kernel whose source does not depend on the PIV geometry. Once per device.

  • initialization() does everything geometry-dependent: builds and uploads the PIV grid and parameter blocks, allocates all device arrays, compiles the two remaining kernels (§5), and plans the batched FFTs. Call it again whenever the image size or the PIV parameters change.

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:

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

select_Platform()

choose the OpenCL platform/device, return the reikna thread

compile_Kernels()

compile all geometry-independent kernels (once per device)

initialization()

grid + buffers + geometry-dependent kernels + FFT plans (once per configuration)

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,” SoftwareX 20, 101256, doi: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,” Exp. Fluids 35, 408–421. — The underlying algorithm (see the PIV algorithm page).

  3. Khronos OpenCL specification — the portable compute standard this implementation targets through pyopencl.