ARKODE: a flexible IVP solver infrastructure for one-step methods

TL;DR

ARKODE introduces a flexible time-integration framework for IVPs of the form $M(t)\,y'(t)=f(t,y)$, supporting explicit RK, diagonally implicit RK, ImEx additive RK, and multirate infinitesimal (MRI) methods. It provides three time-stepping modules—ERKStep, ARKStep, and MRIStep—paired with shared infrastructure for error control, interpolation, event detection, and inequality constraints, enabling rapid experimentation with new methods. Numerical results span simple advection–diffusion–reaction tests to large-scale 3D multiphysics simulations, illustrating performance tradeoffs across solver families (ERKStep, ARKStep in ImEx/DIRK modes, and MRIStep with MRI-GARK/MIS) and highlighting scalable HPC deployments. As part of SUNDIALS, ARKODE offers a robust, extensible platform for developing and deploying high-order, flexible time integrators in multiphysics codes.

Abstract

We describe the ARKODE library of one-step time integration methods for ordinary differential equation (ODE) initial-value problems (IVPs). In addition to providing standard explicit and diagonally implicit Runge--Kutta methods, ARKODE also supports one-step methods designed to treat additive splittings of the IVP, including implicit-explicit (ImEx) additive Runge--Kutta methods and multirate infinitesimal (MRI) methods. We present the role of ARKODE within the SUNDIALS suite of time integration and nonlinear solver libraries, the core ARKODE infrastructure for utilities common to large classes of one-step methods, as well as its use of ``time stepper'' modules enabling easy incorporation of novel algorithms into the library. Numerical results show example problems of increasing complexity, highlighting the algorithmic flexibility afforded through this infrastructure, and include a larger multiphysics application leveraging multiple algorithmic features from ARKODE and SUNDIALS.

Figures (3)

• Figure 1: The relationship between ARKODE and its time-stepping modules. ARKODE provides the time integration loop; each stepper performs individual time steps leveraging the shared infrastructure as needed.
• Figure 2: Work-precision plot for the default 2nd to 5th order ERK methods with various error controllers and tolerances. Configurations further to the left and lower in the plot indicate methods with less work and greater accuracy, respectively.
• Figure 3: Demonstration code scaling: simulation times (first row) and algorithmic scalability (second row). fAdv corresponds to evaluation of the advection operator and fChem the chemistry operator. In the first row, MPI corresponds with the global synchronization times required following each rank-local solver component, and other contains all components that are not directly measured (e.g., time spent within the GPU-based MAGMA linear solver). In the second row, steps, slow steps and fast steps show the total number of time steps taken at each level of each solver, LSetups shows the total number of linear solver setup calls, Jevals shows the total number of Jacobian evaluations, and Niters shows the total number of Newton iterations.