GPSODE: A graphics processing unit-native extrapolation-based implicit stiff ordinary differential equation solver for reacting flow simulations
Description
Accurate simulations of combustion phenomena require solving stiff ordinary differential equations (ODEs) governing chemical reactions, which often constitute a substantial fraction of the total simulation time. However, realistic chemical reaction mechanisms typically involve numerous species and exhibit high levels of mathematical stiffness. While implicit solvers with adaptive step sizes are essential for handling stiff systems, their implementation on Graphics Processing Units (GPUs) faces significant challenges due to thread divergence arising from adaptive time-stepping and error control mechanisms. As a result, existing GPU-accelerated chemical ODE solvers have predominantly focused on explicit methods or simple implicit methods, leaving a significant gap in efficient GPU-accelerated implicit solvers for stiff ODEs requiring high-order accuracy. Here we introduce GPSODE (GPU-accelerated Parallel Stiff ODE solver), a CUDA-centric package (kernels and device-side scheduling in CUDA; Python/C utilities in the repository) that implements an enhanced state-detection-based SEULEX (SD-SEULEX) algorithm specifically optimized for GPU architectures. The SD-SEULEX method enables the software to dynamically adjust the set of active grid points to mitigate GPU thread divergence while maintaining the numerical accuracy and stability of the classical SEULEX method. Floating-point work, state tagging, and index compaction for clustering run on the GPU; the host launches kernels and reads compact counters. Validation tests demonstrate excellent agreement with reference solutions for zero-dimensional homogeneous reactors using H2, CH4, and C7H16 mechanisms, as well as one-dimensional hydrogen-air premixed flames. Performance evaluations show that SD-SEULEX achieves a two-fold performance improvement compared to the baseline SEULEX algorithm. When integrated into OpenFOAM, GPSODE achieves a 172x speedup in the chemical solver compared to a single CPU core, with multi-GPU scaling tests showing 18.4x speedup in the chemical solver and 10.1x overall acceleration versus a 32-core CPU node. C-linkable entry points are provided for integration with OpenFOAM and other CFD stacks.