GPSODE: A graphics processing unit-native extrapolation-based implicit stiff ordinary differential equation solver for reacting flow simulations

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.精准模拟燃烧现象，需求解支配化学反应的刚性常微分方程（Ordinary Differential Equations，简称ODEs），此类方程往往占据总模拟时长的显著比例。然而，实际的化学反应机理通常涉及大量组分，且呈现出极强的数学刚性。采用自适应步长的隐式求解器是处理刚性系统的核心手段，但在图形处理器（Graphics Processing Units，简称GPUs）上部署此类求解器却面临诸多严峻挑战，其根源在于自适应时间步长与误差控制机制引发的线程发散问题。因此，现有的GPU加速化学ODE求解器大多仅支持显式方法或简单隐式方法，在面向高精度刚性ODE的高效GPU加速隐式求解器领域仍存在显著空白。本文提出了GPSODE（GPU加速并行刚性ODE求解器，GPU-accelerated Parallel Stiff ODE solver）：这是一款以CUDA为核心的软件包（代码仓库中包含CUDA内核与设备端调度逻辑，以及Python/C工具集），实现了专为GPU架构优化的、基于状态检测的改进型SEULEX（SD-SEULEX）算法。SD-SEULEX算法允许软件动态调整活跃网格点集合，以缓解GPU线程发散问题，同时保留经典SEULEX算法的数值精度与稳定性。浮点运算、状态标记与聚类索引压缩均在GPU上执行；主机端负责启动内核并读取压缩后的计数器数据。验证测试结果显示，针对采用H₂、CH₄及C₇H₁₆反应机理的零维均质反应器，以及一维氢-空气预混火焰，该求解器的计算结果与参考解吻合极佳。性能评估结果表明，相较于基准SEULEX算法，SD-SEULEX的性能提升可达两倍。将GPSODE集成至OpenFOAM后，相较于单CPU核心，其化学求解器的加速比可达172倍；多GPU扩展性测试显示，化学求解器的加速比为18.4倍，相较于32核CPU节点，整体加速比可达10.1倍。软件提供了可与C语言兼容的接口，便于集成至OpenFOAM及其他计算流体动力学（Computational Fluid Dynamics，简称CFD）框架中。