Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively Refined Meshes for Compressible Fluid Dynamics

TL;DR

The paper develops and analyzes Paired-Explicit Runge-Kutta (P-ERK) methods for multirate time integration on adaptively refined meshes in compressible fluid dynamics. By constructing and optimizing stability polynomials and estimating Jacobian spectra, the authors implement partitioned time stepping on AMR within DGSEM on Trixi.jl, achieving significant speedups on viscous and inviscid problems, with linear stability constraints guiding timestep choices. They provide extensive nonlinear stability analyses, demonstrate performance across a suite of hyperbolic and hyperbolic-parabolic tests, and quantify overheads in computation and memory, along with reproducibility resources. The work advances practical, high-performance time integration for AMR-based simulations, guiding the design of robust multirate schemes for complex flows.

Abstract

In this paper, we apply the Paired-Explicit Runge-Kutta (P-ERK) schemes by Vermeire et. al. (2019, 2022) to dynamically partitioned systems arising from adaptive mesh refinement. The P-ERK schemes enable multirate time-integration with no changes in the spatial discretization methodology, making them readily implementable in existing codes that employ a method-of-lines approach. We show that speedup compared to a range of state of the art Runge-Kutta methods can be realized, despite additional overhead due to the dynamic re-assignment of flagging variables and restricting nonlinear stability properties. The effectiveness of the approach is demonstrated for a range of simulation setups for viscous and inviscid convection-dominated compressible flows for which we provide a reproducibility repository. In addition, we perform a thorough investigation of the nonlinear stability properties of the Paired-Explicit Runge-Kutta schemes regarding limitations due to the violation of monotonicity properties of the underlying spatial discretization. Furthermore, we present a novel approach for estimating the relevant eigenvalues of large Jacobians required for the optimization of stability polynomials.

Figures (13)

• Figure 1: Illustration of the assignment of grid cells (grey boundaries) to partitions $r= 1, 2, 3$ based on the minimum edge length $h$ of the grid cells. Given a base number of stage-evaluations $E^{(1)}$ that is used to integrate the coarse blue cells, a $2E^{(1)} = E^{(2)}$ and a $4E^{(1)} = E^{(3)}$ stage-evaluation method for the green (medium) and orange (fine) cells, respectively, are required to perform the time-integration without reduction in $\text{CFL}$ number (assuming constant wave speeds $\rho_i$ across the domain). The Legendre-Gauss-Lobatto integration points of a second-order ($k=2$) DG method are indicated by the black dots. Note that on each edge there are actually two points stacked on each other, one for each adjacent cell to allow for discontinuous solutions.
• Figure 2: Spectrum $\boldsymbol \sigma(D)$ of the fully discrete matrix $D$ of the $\text{P-ERK}_{2; \{8, 16\}}$ scheme applied to \ref{['eq:UpwindFullSystemODE']}. The unstable region with $\vert \lambda \vert > 1$ is shaded in red. Note that the purely real eigenvalue around $(1, 0)$ is indeed stable with numeric value $1.0 - 28 \varepsilon$ and machine epsilon $\varepsilon = \mathcal{O}\left(10^{-16}\right)$.
• Figure 3: Cell averages $U_i(\Delta t)$ (\ref{['fig:u816_64']}) after performing one timestep of the second-order 8-16 optimized PERK scheme applied to \ref{['eq:UpwindFullSystemODE']} on a uniform mesh with $\Delta x = 2/64$. The entries of the fully discrete Matrix $D$ (cf \ref{['eq:PartitionedSystem']}) are shown in \ref{['fig:MatrixPlot']}.
• Figure 4: Result $\boldsymbol U^1$ to \ref{['eq:PartitionedSystem']} after one timestep of two different 6-12 PERK schemes with different stage evaluation patterns on a uniform mesh. The solution colored in blue correspond to the standard PERK scheme \ref{['eq:PERK_ButcherTableauClassic']} while the solution obtained with the alternating PERK scheme \ref{['eq:PERK_ButcherTableauAlternating']} is colored in orange.
• Figure 5: Reduced spectrum obtained from Arnoldi iterations around shifts (\ref{['fig:ReducedSpectrum']}) and full spectrum obtained from standard eigendecomposition (\ref{['fig:FullSpectrum']}). The spectrum corresponds to a dgsem discretization of the compressible Euler equations with an initial condition leading to a Kelvin-Helmholtz instability.