Seunghyun Joo

Seunghyun Joo

Research Associate at Seoul National University | Looking for a PhD position

CuAD: Source-to-Source Automatic Differentiation for CUDA

GitHub arXiv

🔍 Project Overview

CuAD (CUDA Automatic Differentiation) reads CUDA kernels and writes their derivatives as new CUDA source: a forward (tangent) companion f_d and a reverse (adjoint) companion f_b for each __global__ / __device__ function. Because it differentiates the kernel source itself, every derivative is exact (no finite differences), and the generated kernels run on the GPU like the originals.

CuAD was developed for the exact discrete adjoint of CuDAFLOW, where it differentiates the flux, viscous, and turbulence-model residual kernels.

đź“„ Paper

CuAD: Source-to-Source Automatic Differentiation of CUDA Kernels for Discrete Adjoints in Computational Fluid Dynamics, arXiv:2610.26694 (publication page).

  • Under the no-aliasing convention of CUDA kernels, the reverse sweep recomputes instead of storing. On a Biot–Savart summation the CuAD adjoint costs about four times the primal at every size, while the heap-based tapes of Clad and Enzyme (without __restrict__) cost 9 to 4400 times the primal.
  • For a transonic NACA 0012, the residual, both boundary conditions and the wall force are differentiated without change to the source, and the adjoint and tangent shape sensitivities agree to 1e-13.
Biot-Savart adjoint cost of CuAD, Clad and Enzyme
Biot–Savart adjoint of four tools under one compiler: cost relative to the primal kernel and the device heap required by the tape or cache.
Transonic NACA 0012 lift adjoint and shape sensitivities
Transonic NACA 0012: lift adjoint and the lift and drag shape sensitivities from the adjoint (bars) and the tangent (circles).

🛠️ Key Features

  • Forward and reverse mode, including calls between user functions (each callee gets its own companion).
  • Activity analysis: only variables that depend on active inputs are differentiated.
  • Recompute rather than store: under the convention that a kernel’s inputs and outputs do not alias, loops are recomputed in the reverse sweep; only scalars carried between iterations are checkpointed, and pure accumulations (v = v + e) need nothing.
  • Thread-local tape: values that must be restored go to a fixed-size per-thread array, not to the device heap.
  • GPU-safe adjoint accumulation: atomicAdd for shared arrays, plain += for thread-local data.
  • Loops and control flow: grid-stride loops, per-iteration recomputation, #pragma unroll, early exits, c ? a : b.
  • C++ constructs: function templates, default arguments, reference outputs, local pointers and arrays of pointers.
  • Math library: trigonometric, exponential, power, atan2, hypot, copysign, cbrt, and more, in float and double.

🆕 Improvements over the Version Bundled in CuDAFLOW

  • Fixes for silently wrong derivatives, each checked against finite differences: c ? a : b (its adjoint was dropped), a one-sided if in the middle of a function (the reverse sweep ran in the wrong order), value-returning user functions nested in expressions (adjoint was 0), and copysign.
  • Less tape: pure accumulations need no tape, and recompute loops replay their primal without pushes.
  • Broader language coverage: templates, default arguments, reference outputs, array parameters, local pointers, early exits, branch-local declarations, bare blocks, #pragma unroll, in-place updates of local arrays in loops.
  • More math functions: atan2, atan, asin, acos, hypot, cbrt, copysign, and piecewise-constant floor / ceil / round.

âś… Verification

The test suite checks every generated tangent and adjoint against central finite differences and with adjoint dot-product tests (⟨y̅, J ẋ⟩ = ⟨Jᵀ y̅, ẋ⟩):

CaseCheck
Small kernels (scalar math, loops, overwritten values, device calls, atomic scatter)Tangent and adjoint vs. FD, relative error ~1e-11
CuDAFLOW kernels (conservative-to-primitive, cell volume)Tangent FD and adjoint dot test
Mesh shape optimisation and NACA 0012 inverse designAdjoint gradient vs. FD, then converged optimisation
Spalart–Allmaras field inversion (FIML)SA source adjoint dot test; recovers the correction field β
Loops carrying scalars between iterations (v0.2.1)Adjoint vs. FD, relative error ~1e-10
Transonic NACA 0012, 2-D Euler (paper)Adjoint vs. tangent shape sensitivities, 5e-14

🧪 Technical Skills Demonstrated

  • Languages / tools: C++17, CUDA, Clang LibTooling (AST analysis and source rewriting)
  • Methods: Algorithmic differentiation (tangent and adjoint), activity analysis, taping, discrete adjoint methods

GPL-3.0.