Additive Runge-Kutta Research Articles

A performance study of horizontally explicit vertically implicit (HEVI) time-integrators for non-hydrostatic atmospheric models

We conduct a thorough study of different forms of horizontally explicit and vertically implicit (HEVI) time-integration strategies for the compressible Euler equations on spherical domains typical of nonhydrostatic global atmospheric applications. We compare the computational time and complexity of two nonlinear variants (NHEVI-GMRES and NHEVI-LU) and a linear variant (LHEVI). We report on the performance of these three variants for a number of additive Runge-Kutta methods ranging in order of accuracy from second through fifth, and confirm the expected order of accuracy of the HEVI methods for each time-integrator. To gauge the maximum usable time-step of each HEVI method, we run simulations of a nonhydrostatic baroclinic instability for 100 days and then use this time-step to compare the time-to-solution of each method. The results show that NHEVI-LU is 3x faster than NHEVI-GMRES, and LHEVI is 8x faster than NHEVI-GMRES, for the idealized cases tested. The baroclinic instability and inertia-gravity wave simulations indicate that the optimal choice of HEVI method is LHEVI. For the fastest time-to-solution, the second and third order methods are best although the better accuracy of the high-order methods (particularly the fourth order method) should be considered.In the future, we will report on whether these results hold for more complex problems using, e.g., real atmospheric data and/or a higher model top typical of space weather applications.

Implicit–explicit multirate infinitesimal stage-restart methods

Implicit–Explicit (IMEX) methods are flexible numerical time integration methods which solve an initial-value problem (IVP) that is split into stiff and nonstiff processes with the goal of lower computational costs than a purely implicit or explicit approach. A complementary form of flexible IVP solvers are multirate infinitesimal methods for problems split into fast- and slow-changing dynamics, that solve a multirate IVP by evolving a sequence of “fast” IVPs using any suitably accurate algorithm. This article introduces a new class of high-order implicit–explicit multirate methods that are designed for multirate IVPs in which the slow-changing dynamics are further split in an IMEX fashion. This new class, which we call implicit–explicit multirate infinitesimal stage-restart (IMEX-MRI-SR), both improves upon the previous implicit–explicit multirate infinitesimal generalized-structure additive Runge Kutta (IMEX-MRI-GARK) methods by allowing for far easier creation of new embedded methods, and extends multirate exponential Runge Kutta (MERK) methods by allowing the fast-changing dynamics to be nonlinear and the methods to be implicit. We leverage GARK theory to derive conditions for orders of accuracy up to four, and we provide second- and third-order accurate example methods, which are the first known embedded MRI methods with IMEX structure. We then perform numerical simulations demonstrating convergence rates and computational performance in both fixed-step and adaptive-step settings.

ARKODE: A Flexible IVP Solver Infrastructure for One-step Methods

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 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.

An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement

We present an accurate and efficient solver for atmospheric dynamics simulations that allows for non-conforming mesh refinement. The model equations are the conservative Euler equations for compressible flows. The numerical method is based on an h−adaptive Discontinuous Galerkin spatial discretization and on a second order Additive Runge Kutta IMEX method for time discretization, especially designed for low Mach regimes. The solver is implemented in the framework of the deal.II library, whose mesh refinement capabilities are employed to enhance efficiency. A number of numerical experiments based on classical benchmarks for atmosphere dynamics demonstrate the properties and advantages of the proposed method.

Fractional-step Runge–Kutta methods: Representation and linear stability analysis

Fractional-step methods are a popular and powerful divide-and-conquer approach for the numerical solution of differential equations. When the integrators of the fractional steps are Runge–Kutta methods, such methods can be written as generalized additive Runge–Kutta (GARK) methods, and thus the representation and analysis of such methods can be done through the GARK framework. We show how the general Butcher tableau representation and linear stability of such methods are related to the coefficients of the splitting method, the individual sub-integrators, and the order in which they are applied. We use this framework to explain some observations in the literature about fractional-step methods such as the choice of sub-integrators, the order in which they are applied, and the role played by negative splitting coefficients in the stability of the method.

An efficient IMEX-DG solver for the compressible Navier-Stokes equations for non-ideal gases

We propose an efficient, accurate and robust IMEX solver for the compressible Navier-Stokes equations describing non-ideal gases with general cubic equation of state and Stiffened-Gas EOS. The method is based on an h-adaptive Discontinuous Galerkin spatial discretization and on an Additive Runge Kutta IMEX method for time discretization. It is specifically tailored for low Mach number applications and allows to simulate low Mach regimes at a significantly reduced computational cost, while maintaining full second order accuracy also for higher Mach number regimes. The method has been implemented in the framework of the deal.II numerical library, whose adaptive mesh refinement capabilities are employed to enhance efficiency. Refinement indicators appropriate for real gas phenomena have been introduced. A number of numerical experiments on classical benchmarks for compressible flows and their extension to real gases demonstrate the properties of the proposed method.

A Fast Time-Stepping Strategy for Dynamical Systems Equipped with a Surrogate Model

Simulation of complex dynamical systems arising in many applications is computationally challenging due to their size and complexity. Model order reduction, machine learning, and other types of surrogate modeling techniques offer cheaper and simpler ways to describe the dynamics of these systems, but are inexact and introduce additional approximation errors. In order to overcome the computational difficulties of the full complex models, on one hand, and the limitations of surrogate models, on the other, this work proposes a new accelerated time-stepping strategy that combines information from both. This approach is based on the multirate infinitesimal general-structure additive Runge-Kutta (MRI-GARK) framework. The inexpensive surrogate model is integrated with a small timestep to guide the solution trajectory, and the full model is treated with a large timestep to occasionally correct for the surrogate model error and ensure convergence. We provide a theoretical error analysis, and several numerical experiments, to show that this approach can be significantly more efficient than using only the full or only the surrogate model for the integration.

High‐order implicit‐explicit additive Runge–Kutta schemes for numerical combustion with adaptive mesh refinement

AbstractThe objective of this study is to develop and apply efficient solution techniques for numerical modeling of combustion with stiff chemical kinetics in practical combustors. The new technique combines a fourth‐order implicit‐explicit (IMEX) additive Runge‐Kutta scheme (ARK) with adaptive mesh refinement (AMR). The IMEX component treats the stiff reactions implicitly but integrates convection and diffusion explicitly in time, and thus permits the solution to advance with larger time‐step sizes than that of explicit time‐marching methods alone. The AMR further adds computational efficiency by effectively placing high spatial resolution meshes in regions with strong gradients, such as flame fronts. The novelty of this study is in the integration of a fourth‐order IMEX ARK method with AMR for a high‐order finite‐volume scheme and the application to solving complex reacting flows governed by the compressible Navier–Stokes equations with very stiff chemistry in a practical combustor geometry. The effectiveness and performance of the adaptive ARK4 is assessed for complex reacting flows by examining properties, such as the presence of shock waves, the time‐scale changes in response to AMR levels, and the size and stiffness of reaction mechanisms for various fuels such as , , and . The new adaptive ARK4 method is verified and validated using a convection‐diffusion‐reaction problem and shock‐driven combustion, respectively. The validated algorithm is then applied to solve the stiff ‐air combustion in a bluff‐body combustor. A significant speedup of three orders of magnitude is achieved in comparison to the standard ERK4 method at the given solution accuracy.

Symplectic P-stable additive Runge—Kutta methods

<p style='text-indent:20px;'>Classical symplectic partitioned Runge–Kutta methods can be obtained from a variational formulation where all the terms in the discrete Lagrangian are treated with the same quadrature formula. We construct a family of symplectic methods allowing the use of different quadrature formulas (primary and secondary) for different terms of the Lagrangian. In particular, we study a family of methods using Lobatto quadrature (with corresponding Lobatto IIIA-B symplectic pair) as a primary method and Gauss–Legendre quadrature as a secondary method. The methods have the same implicitness as the underlying Lobatto IIIA-B pair, and, in addition, they are <i>P-stable</i>, therefore suitable for application to highly oscillatory problems.</p>

A unified formulation of splitting-based implicit time integration schemes

Splitting-based time integration approaches such as fractional step, alternating direction implicit, operator splitting, and locally one dimensional methods partition the system of interest into components, and solve individual components implicitly in a cost-effective way. This work proposes a unified formulation of splitting time integration schemes in the framework of general-structure additive Runge–Kutta (GARK) methods. Specifically, we develop implicit-implicit (IMIM) GARK schemes, provide the order conditions for this class, and explain their application to partitioned systems of ordinary differential equations. We show that classical splitting methods belong to the IMIM GARK family, and therefore can be studied in this unified framework. New IMIM-GARK splitting methods are developed and tested using parabolic systems.

Semi-Implicit and Explicit Runge Kutta Methods for Stiff Ordinary Differential Equations

In this work, we study the A[α] – stability of the additive methods of Runge- Kutta kind of orders ranging from 2 up to 4 that will be applied for determining some stiff nonlinear system of the ODEs. Moreover, we find the stability function for the additive Runge-Kutta method and some methods of this type of order 2,3, and 4. Where the method (A,B1 ) is A-stable and semi-implicit and method (A,B2 ) is explicit. Furthermore, the stiff term is managed by the semi-implicit Runge-Kutta method while no stiff term is treated by the explicit Runge Kutta method. Those methods are suitable for solving chemical reactions problems that include stiff and non-stiff terms.

Linearly implicit GARK schemes

Systems driven by multiple physical processes are central to many areas of science and engineering. Time discretization of multiphysics systems is challenging, since different processes have different levels of stiffness and characteristic time scales. The multimethod approach discretizes each physical process with an appropriate numerical method; the methods are coupled appropriately such that the overall solution has the desired accuracy and stability properties. The authors developed the general-structure additive Runge–Kutta (GARK) framework, which constructs multimethods based on Runge–Kutta schemes.This paper constructs the new GARK-ROS/GARK-ROW families of multimethods based on linearly implicit Rosenbrock/Rosenbrock-W schemes. For ordinary differential equation models, we develop a general order condition theory for linearly implicit methods with any number of partitions, using exact or approximate Jacobians. We generalize the order condition theory to two-way partitioned index-1 differential-algebraic equations. Applications of the framework include decoupled linearly implicit, linearly implicit/explicit, and linearly implicit/implicit methods. Practical GARK-ROS and GARK-ROW schemes of order up to four are constructed.

Numerical investigation of self-sustaining modes of 2D planar detonations under concentration gradients in hydrogen–oxygen mixtures

Self-sustaining modes of detonations in non-uniform mixtures for hydrogen and oxygen are numerically studied. Concentration gradients are set vertical to the direction of detonation propagation, with equivalence ratios decreasing from rich values to lean ones. High-resolution codes based on the fifth-order weighted essentially non-oscillatory (WENO) scheme in spatial discretization and the 3rd-order Additive Runge-Kutta (ARK) schemes in time discretization, with detailed chemical kinetics model, are developed. Two-dimensional (2D) numerical simulations are conducted in a 2D planar geometry, non-uniform distributions are then precisely controlled as a function of width. The cellular structures and transmissions, failure and re-initiation mechanisms of detonation fronts are characterized. Results show that the low concentration gradient produces the multi-head mode of detonation with continuous adaptation of velocities and multicellular structures to local concentration. The high concentration gradient generates single-head mode with detonation re-initiation at channel walls produced by reflected transverse Mach waves. Specially, the intermediate distribution between the low and high concentration gradients induces an alternation mode between singe-head and multi-head detonations.

Effect of transversal concentration gradient on H2–O2 cellular detonation**Project supported by the Natural National Science Foundation of China (Grant Nos. 11972090, 11732003, and U1830139), the Beijing Natural Science Foundation, China (Grant No. 8182050), and the National Key Research and Development Program of China (Grant No. 2017YFC0804700).

A two-dimensional detonation in H2–O2 system is simulated by a high-resolution code based on the fifth-order weighted essentially non-oscillatory (WENO) scheme in the spatial discretization and the 3th-order additive Runge–Kutta schemes in the time discretization, by using a detailed chemical model. The effect of a concentration gradient on cellular detonation is investigated. The results show that with the increase of the concentration gradient, the cell instability of detonation increases and gives rise to the oscillation of average detonation velocity. After a long time, for the case of the lower gradient the detonation can be sustained, with the multi-head mode and single-head mode alternating, while for the higher gradient it propagates with a single-head mode.

Coupled Multirate Infinitesimal GARK Schemes for Stiff Systems with Multiple Time Scales

Traditional time discretization methods use a single timestep for the entire system of interest and can perform poorly when the dynamics of the system exhibits a wide range of time scales. Multirate infinitesimal step (MIS) methods (Knoth and Wolke, 1998) offer an elegant and flexible approach to efficiently integrate such systems. The slow components are discretized by a Runge-Kutta method, and the fast components are resolved by solving modified fast differential equations. Sandu (2018) developed the Multirate Infinitesimal General-structure Additive Runge-Kutta (MRI-GARK) family of methods that includes traditional MIS schemes as a subset. The MRI-GARK framework allowed the construction of the first fourth order MIS schemes. This framework also enabled the introduction of implicit methods, which are decoupled in the sense that any implicitness lies entirely within the fast or slow integrations. It was shown by Sandu that the stability of decoupled implicit MRI-GARK methods has limitations when both the fast and slow components are stiff and interact strongly. This work extends the MRI-GARK framework by introducing coupled implicit methods to solve stiff multiscale systems. The coupled approach has the potential to considerably improve the overall stability of the scheme, at the price of requiring implicit stage calculations over the entire system. Two coupling strategies are considered. The first computes coupled Runge-Kutta stages before solving a single differential equation to refine the fast solution. The second alternates between computing coupled Runge-Kutta stages and solving fast differential equations. We derive order conditions and perform the stability analysis for both strategies. The new coupled methods offer improved stability compared to the decoupled MRI-GARK schemes. The theoretical properties of the new methods are validated with numerical experiments.

The Bloch equation for spin dynamics in electron storage rings: Computational and theoretical aspects

In this paper, we describe our work on spin polarization in high-energy electron storage rings which we base on the Full Bloch equation (FBE) for the polarization density and which aims towards the [Formula: see text] option of the proposed Future Circular Collider (FCC-ee) and the proposed Circular Electron Positron Collider (CEPC). The FBE takes into account non spin-flip and spin-flip effects due to synchrotron radiation including the spin-diffusion effects and the Sokolov–Ternov effect with its Baier–Katkov generalization as well as the kinetic-polarization effect. This mathematical model is an alternative to the standard mathematical model based on the Derbenev–Kondratenko formulas. For our numerical and analytical studies of the FBE, we develop an approximation to the latter to obtain an effective FBE. This is accomplished by finding a third mathematical model based on a system of stochastic differential equations (SDEs) underlying the FBE and by approximating that system via the method of averaging from perturbative ODE theory. We also give an overview of our algorithm for numerically integrating the effective FBE. This discretizes the phase space using spectral methods and discretizes time via the additive Runge–Kutta (ARK) method which is a high-order semi-implicit method. We also discuss the relevance of the third mathematical model for spin tracking.

Extended multirate infinitesimal step methods: Derivation of order conditions

Multirate methods are specially designed for problems with multiple time scales. The multirate infinitesimal step method (MIS) was developed as a generalization of the so called split-explicit Runge–Kutta methods, where the integration of the fast part is conducted analytically. The MIS method was originally evolved for applications related to numerical weather prediction, i.e. the integration of the compressible Euler equation.In this work, an extension to MIS methods will be presented where an arbitrary Runge–Kutta method (RK) is applied for the integration of the fast component. Furthermore, the order convergence from the original MIS method will be reinvestigated including the derivation of conditions up to order four. Additionally will be presented how well-known methods such as recursive flux splitting multirate method, (Schlegel et al., 2012) partitioned Runge–Kutta method, (Jackiewicz and Vermiglio, 2000) and generalized additive Runge–Kutta method, (Sandu and Günther, 2015) are related to or can be cast as an extended MIS method. An exemplary MIS method of order four with five stages will show that the convergence behaviour not only depends on the applied method for the integration of the fast component. The method will further indicate that the used fast time step plays a significant role.

Diagonally implicit Runge–Kutta methods for stiff ODEs

Based principally on a recent review of diagonally implicit Runge–Kutta (DIRK) methods applied to stiff first-order ordinary differential equations (ODEs) by the present authors, several nearly optimal, general purpose, DIRK-type methods are presented. Methods given range from third- to sixth-order in four- to nine-stages. All are both stiffly-accurate and L-stable and many are internally L-stable on stages where this is possible. Focus is placed on stage-order two methods. To facilitate step-size control via error estimation, an embedded method is included with each method listed. As the present methods must make no concessions to an explicit method as part of an additive Runge–Kutta method, the present ESDIRKs are generally better than their IMEX analogs. Five new explicit first-stage, singly diagonally-implicit Runge–Kutta (ESDIRK) methods are presented based on lessons learned from the review. To the best of our knowledge, the sixth-order, L-stable, stage-order two, 6(5)-pair, ESDIRK included herein is the only one of its kind.

A Class of Multirate Infinitesimal GARK Methods

Differential equations arising in many practical applications are characterized by multiple time scales. Multirate time integration seeks to solve them efficiently by discretizing each scale with a different, appropriate time step, while ensuring the overall accuracy and stability of the numerical solution. In a seminal paper Knoth and Wolke (APNUM, 1998) proposed a hybrid solution approach: discretize the slow component with an explicit Runge-Kutta method, and advance the fast component via a modified fast differential equation. The idea led to the development of multirate infinitesimal step (MIS) methods by Wensch et al. (BIT, 2009.)G\"{u}nther and Sandu (BIT, 2016) explained MIS schemes as a particular case of multirate General-structure Additive Runge-Kutta (MR-GARK) methods. The hybrid approach offers extreme flexibility in the choice of the numerical solution process for the fast component. This work constructs a family of multirate infinitesimal GARK schemes (MRI-GARK) that extends the hybrid dynamics approachin multiple ways. Order conditions theory and stability analyses are developed, and practical explicit and implicit methods of up to order four are constructed. Numerical results confirm the theoretical findings. We expect the new MRI-GARK family to be most useful for systems of equations with widely disparate time scales, where the fast process is dispersive, and where the influence of the fast component on the slow dynamics is weak.

Higher-order additive Runge–Kutta schemes for ordinary differential equations

Two new implicit–explicit, additive Runge–Kutta (ARK2) methods are given with fourth- and fifth-order formal accuracies, respectively. Both combine explicit Runge–Kutta (ERK) methods with explicit, singly-diagonally implicit Runge–Kutta (ESDIRK) methods and include an embedded method for error control. The two methods have ESDIRKs which are internally L-stable on stages three and higher, have only modestly negative eigenvalues to the stage and step algebraic-stability matrices and have stage-order two. To improve computational efficiency, the fourth-order method has a diagonal coefficient of 0.1235. This is done to offset much of the extra computational cost of an extra stage by facilitating iterative convergence at each stage. Linear stability domains for both ERK methods have been made quite large and the dominant coupling stability term between the stability of the ESDIRK and ERK for very stiff modes has been removed. Though the fourth-order method is one of the best all-around fourth-order IMEX ARK2 of which we are aware, the fifth-order method is likely best suited to mildly stiff problems with tight error tolerances. Methods are tested using the Van der Pol and Kaps' singular-perturbation problems. Results suggest that these new methods represent an improvement over existing methods of the same class.