New Fourth-Order Time-Stepping Schemes Promise Faster Flow Simulations

Numerical simulations of flows, sound waves, and other time-dependent physical phenomena often grind along at the pace of their slowest region. A team of applied mathematicians and computational engineers now reports a way to let different parts of a simulation advance at different speeds without sacrificing accuracy, and the result is a family of fourth-order accurate time-integration schemes that can cut computing time dramatically in the right circumstances. The work, published in the open-access journal Computational Science and Engineering, extends a class of methods known as Paired-Explicit Runge-Kutta (P-ERK) schemes to fourth order of consistency, a level of accuracy that had previously been out of reach for this approach.

The researchers, Daniel Doehring, Lars Christmann, Michael Schlottke-Lakemper, Gregor J. Gassner, and Manuel Torrilhon, built their new schemes on the framework of partitioned Runge-Kutta methods, which apply different integration rules to different components of a system of equations. The original P-ERK schemes were introduced by Vermeire and colleagues in 2019 and extended to third order in 2022. The new contribution pushes that framework to fourth order, which required overcoming a substantial mathematical obstacle: the order conditions for fourth-order partitioned methods include coupling constraints that tie the coefficient tables of all the individual schemes together, making their construction significantly more difficult than at lower orders.

The core problem the methods address is the Courant-Friedrichs-Lewy (CFL) condition, a fundamental restriction on explicit time stepping. When a simulation uses an explicit method, the size of each time step is limited by the fastest wave speed anywhere in the computational domain. On a non-uniform grid, where some cells are tiny and others are large, or in a problem with multiscale physics, the tiniest cells or fastest regions dictate the step size for the entire simulation. Large cells that could safely take much bigger steps are forced to march in lockstep with the smallest ones, wasting computational effort on work that is far more accurate than it needs to be.

Multirate and local time-stepping methods have been developed since the 1980s to relax this constraint by allowing different timesteps in different regions. The P-ERK approach takes a particularly elegant route: it pairs together several optimized explicit Runge-Kutta schemes, each with a different stability domain, and applies them to different partitions of the computational domain. Crucially, the schemes can be implemented much like well-known implicit-explicit methods, requiring no changes to the spatial discretization and no specialized interpolation of interface values. The authors emphasize that their schemes work seamlessly with codes employing a method-of-lines approach, in which space and time are discretized separately.

The path to fourth order was blocked by the coupling conditions among the coefficient tables, known as Butcher arrays, of the different schemes. Rather than optimizing an entire family of arrays simultaneously, which would create a vastly more complex optimization problem, the authors chose a specific functional form of the arrays that satisfies the coupling conditions by construction. This choice allows each method in the family to be optimized separately and makes it easy to extend the family with additional schemes. The resulting structure leads to a special form of the stability polynomials, the mathematical objects that govern how the method damps numerical errors, which must be respected during optimization of the stability domain.

A notable consequence of the design is that a fourth-order P-ERK method requires at least five stages, whereas second- and third-order versions could be built with two or three stages. The authors also quantify the price of their design choices: the constrained stability polynomials yield smaller maximum admissible timesteps than completely free ones, but the loss is modest. For a six-stage polynomial with one free coefficient, roughly 85 percent of the unconstrained timestep can be achieved, and this ratio rises to about 96.5 percent for a sixteen-stage polynomial. The team validated the new schemes extensively, demonstrating linear stability, internal consistency, the designed fourth-order convergence, and conservation of linear invariants such as mass, momentum, and energy.

The practical payoff appears in a series of demanding applications. The authors applied the multirate schemes to inviscid and viscous problems with locally varying wave speeds, induced either by non-uniform grids or by multiscale properties of the governing partial differential equations. Test cases included the linearized Euler equations in three dimensions with a non-uniform mean speed of sound, isentropic vortex advection on adaptively refined meshes, laminar flow around an SD7003 airfoil, and the sound generation of two spinning, co-rotating vortices in a coupled flow-acoustics simulation. Compared to state-of-the-art optimized standalone methods, the multirate P-ERK schemes allowed significant reductions in right-hand-side evaluations and wall-clock time, ranging from 66 percent up to factors greater than four.

One of the most striking demonstrations involved a three-dimensional linearized Euler simulation with a spatially varying speed of sound, discretized with more than 100 million degrees of freedom. The authors constructed ten P-ERK schemes with different numbers of stage evaluations, matched to the local speed of sound, and compared them against several established fourth-order integrators. The multirate scheme outperformed them all in runtime and in the number of scalar right-hand-side evaluations. The authors also note that savings in function evaluations do not translate exactly proportionally into runtime savings, because their solver is memory-bound and the partitioned approach loses some data locality, an honest accounting that adds credibility to the reported speedups.

The schemes also proved well suited to parallel computing. The authors demonstrated speedups for single-threaded program executions, shared-memory parallelism with multiple threads, and distributed-memory parallelism using MPI. For distributed-memory runs, the partitioning of unknowns across ranks can be done independently on each rank, and the authors describe a load-balancing strategy that balances the number of cells per integration level among ranks, which improved performance for the multirate schemes. A reproducibility repository is provided that enables readers to examine all results presented in the work, and the authors indicate that future work will explore P-ERK schemes with error-based timestep control.

For the computational science community, the significance of this work lies in combining three qualities that are often in tension: high order of accuracy, multirate flexibility, and ease of implementation within existing method-of-lines codes. The fourth-order P-ERK schemes satisfy linear stability, internal consistency, designed order of convergence, and conservation of linear invariants while remaining non-intrusive, requiring no modification of the spatial discretization.

Subject of Research: Fourth-order multirate explicit Runge-Kutta time integration methods for partial differential equations

Article Title: Fourth-order paired-explicit Runge-Kutta methods

Article References: Doehring, D., Christmann, L., Schlottke-Lakemper, M., Gassner, G. J., & Torrilhon, M. (2025). Fourth-order paired-explicit Runge-Kutta methods. Computational Science and Engineering, 2(1), Article 1. https://doi.org/10.1007/s44207-025-00005-4

Image Credits: AI Generated

DOI: 10.1007/s44207-025-00005-4

Keywords: numerical analysis, Runge-Kutta methods, multirate time integration, partial differential equations, computational fluid dynamics, CFL condition, stability polynomials, discontinuous Galerkin, scientific computing, high-order methods, parallel computing, flow simulation