CUDA (GPU) backend

Header: cuda_solver.cuh. Namespace: adrt::cuda. Enable with -DADRT_ENABLE_CUDA=ON at configure time.

The CUDA backend is a batched solver: it processes an entire spectrum in parallel with one thread per wavenumber. It provides template-specialised kernels (solveKernel<N>) for \(N = 2, 4, 8, 16\) with register-resident matrix operations, and a cuBLAS-based path for larger \(N\). Linear systems are factorised in double precision.

Configuration

struct adrt::cuda::BatchConfig

Scalars shared across all wavenumbers. Mirrors ADConfig but with the batch dimension num_wavenumbers and a num_moments_max (the maximum number of Legendre moments across layers). Fields include num_layers, num_quadrature (2, 4, 8, 16, or 32), use_delta_m, use_thermal_emission, use_diffusion_lower_bc, surface_albedo, surface_temperature, top_temperature, solar_flux, solar_mu, wavenumber_low / wavenumber_high, and

double spectrum_scaling: A constant factor applied to flux_up only — for example \(10^{-3}(R_p/R_s)^2\) when producing a planet/star flux ratio.

Device data layouts

struct adrt::cuda::DeviceData: Caller-owned device pointers using a structure-of-arrays layout, array[wav * stride + layer_or_level]. Inputs include delta_tau, single_scat_albedo, phase_moments (with phase_moments_shared if the same moments apply to every wavenumber), and either temperature (for on-device Planck evaluation) or planck_levels. Outputs are flux_up, flux_down, and optional flux_direct.

struct adrt::cuda::RawDeviceData: A raw-input interface for retrieval codes such as BeAR. Takes level-major coefficient arrays with BOA at index 0 (absorption_coeff, scattering_coeff, altitude, temperature, wavenumber, and an optional cloud_optical_depth). The solver computes optical depths (trapezoidal rule), single-scattering albedos, and the single-wavenumber Planck function inline, and reverses the ordering to TOA-first internally.

Entry points

void adrt::cuda::solveBatch(const BatchConfig &config, const DeviceData &data, cudaStream_t stream = 0): Launch the batched kernel; all data must already be on the GPU. Call uploadQuadratureData() once beforehand.

HostResult adrt::cuda::solveBatchHost(const BatchConfig &config, const std::vector<float> &delta_tau, const std::vector<float> &single_scat_albedo, const std::vector<float> &phase_moments, bool phase_moments_shared, const std::vector<float> &planck_levels, const std::vector<float> &temperature = {}): Convenience wrapper for standalone use: allocates device memory, copies the host arrays over, runs the kernel, and copies the results back into a HostResult (flux_up, flux_down, flux_direct).

void adrt::cuda::solveBatchFromCoefficients(const BatchConfig &config, const RawDeviceData &data, cudaStream_t stream = 0): Solve from raw retrieval inputs for \(N \leq 8\) with zero intermediate allocation — optical properties and Planck values are computed inside the kernel. Suitable for very large wavenumber counts.

void adrt::cuda::solveBatchFromCoefficients(const BatchConfig &config, const RawDeviceData &data, SolverWorkspaceGPU &workspace, cudaStream_t stream = 0): The \(N > 8\) overload: preprocesses the raw inputs into intermediate arrays held by a persistent SolverWorkspaceGPU, then runs the cuBLAS-based solver.

struct adrt::cuda::SolverWorkspaceGPU: Persistent GPU workspace for the cuBLAS path (\(N > 8\)). Allocate once with allocate(nwav, nlay) and reuse across retrieval iterations; it frees its device memory on destruction.

Note

The CUDA backend uses single precision for the stored optical properties and double precision for the LU factorisation of the linear systems. Analytic temperature Jacobians are not available on this backend.