A new domain‐based implicit‐explicit time stepping scheme based on the class of exponential integrators called sEPIRK

Reliable simulations of flows in real applications involve the task of discretizing both space and time in an accurate and efficient way. To cope with the large semidiscrete systems resulting from a space discretization on appropriate grids, which often include comparatively few very small cells near solid walls for boundary layer resolution, an implicit‐explicit time stepping scheme can be the most efficient variant. We present such a type of scheme utilizing recent exponential integrators called sEPIRK and show its computational advantages opposed to IMEX‐Runge‐Kutta schemes.


Introduction
The spatial discretization of PDEs such as the Navier-Stokes-equations results in a large, stiff system of ODEs which is w.l.o.g. autonomous. The objective of the domain-based IMEX setting is the combination of an IMplicit and an EXplicit time stepping scheme to loosen the strict CFL-condition of the explicit method on the one hand, and reducing the sizes of the nonlinear and linear systems of the implicit scheme on the other hand. To this end, the computational domain is subdivided into an implicit region linked to implicit time stepping containing the small cells and an explicit region linked to explicit time stepping containing the remaining larger cells, see [1]. The components connected to the implicit region can be collected in U I and the other components in U E , so that the ODE system can be re-arranged such that F =F E +F I with: .
Kanevsky et. al. [1] showed good performance of IMEX-Runge-Kutta (IMEX-RK) methods in this domain-based IMEX setting. Our goal is to adopt the exponential integrators of Tokman into this domain-based IMEX-setting, since they turned out to be perfectly suited for the application to large, stiff ODE systems instead of conventional implicit integrators, see [5].

IMEX-(s)EPIRK schemes
In [3] we demonstrated that the additive coupling of explicit Runge-Kutta methods with EPIRK schemes is limited to first order of convergence and thus constructed the simple first order IMEX-EPIRK1 scheme, which we will call IMEX-EPIRK in the following. Further, to obtain higher order schemes we introduced the idea of utilizing the sEPIRK schemes [4], which are tailored for the special case that the right-hand side F can be split as F(U(t)) = LU(t)+N(U(t)), L ∈ R N ×N . We first divide F according to (2) and then define L := DF I (U n ), U n ≈ U(t n ) and N(U) := F(U) − LU at each time step n. Applying the sEPIRK schemes to this domain-based IMEX right-hand side, our IMEX-sEPIRK schemes result as presented in [3]. We can show that the IMEX-sEPIRK schemes in fact consist of coupling EPIRK schemes for the implicit region with explicit Runge-Kutta methods for the explicit region. Taking the coefficient sets for sEPIRK3 and sEPIRK3b from [4] we arrive at the third order IMEX-sEPIRK3 and IMEX-sEPIRK3b schemes.

Numerical experiments
In this section we want to contrast our first order IMEX-EPIRK to the first order IMEX-RK 1 method, which consists of coupling the explicit and the implicit Euler method, and our third order IMEX-sEPIRK3(b) schemes to the IMEX-RK 3 method derived in [2] under the name ARK3(2)4L [2]SA. The IMEX-sEPIRK schemes employ an efficient adaptive Krylov algorithm for the approximation of the ϕ-times-vector operations, which is available in the EPIC package at http://faculty. ucmerced.edu/mtokman/#software. In our implementation we slightly modified that algorithm to primarily allow for a fair 2 of 2 Section 18: Numerical methods of differential equations comparison with the Newton-routines within the IMEX-RK schemes, i.e. enable conformity of the maximum dimension of the Krylov subspaces, the maximum number of iterations and the tolerance. Moreover, small bugfixes were incorporated.

1D viscous Burgers equation
To begin with, we consider the FD-discretized 1D viscous Burgers equation with a steep tangent wave as initial condition, see [1]. We split the domain Ω = [−1, 1] into an explicit region with cell size ∆x E = 10 −3 and an implicit region located at the wave position with cell size ∆x I = ∆x E /S dependent on the stiffness factor S.
First, we verify the experimental order of convergence of the new IMEX-schemes, see Fig. 1 for S = 6 as an example. Next, CPU-Precision diagrams for various stiffness factors S = 6, 12, 24, 48 with CFL-numbers between 51.2 and 0.0005 have been generated. The results show a comparable performance of the IMEX-SEPIRK3(b) schemes to IMEX-RK 3 for S = 6, but an increasing efficiency gain of the IMEX-sEPIRK3(b) schemes with increasing stiffness factor S. In Fig. 2 we see a factor of up to 10 between the IMEX-sEPIRK3(b) and IMEX-RK 3 and a factor of up to 16 between the IMEX-EPIRK and the IMEX-RK 1 for the same accuracy at S = 48. We will see in the Shu-Vortex test case below, that the explanation for that efficiency gain mainly lies in the higher involvement of the exponential integrators within the IMEX-(s)EPIRK schemes due to larger implicit regions for increasing S (6% implicit cells for S = 6 and approximately 36% for S = 48).
Shu-Vortex test case On a periodic triangulation the 2D Euler equations are discretized with DG-P2 for the Shu-Vortex initial condition. The implicit region was set on the refined trace of the traveling vortex, resulting in approximately 50% of implicit cells. In that case the new IMEX-sEPIRK3(b) schemes were 1.8 times faster than IMEX-RK 3, whereas the IMEX-RK 1 showed better performance than the IMEX-EPIRK with a factor of 5.
Furthermore, we analyze the behaviour of the implicit and the explicit solver within the IMEX-schemes by treating the whole domain implicitly resp. explicitly. It turns out, that the explicit RK methods within the IMEX-sEPIRK3(b) schemes perform comparable to that within the IMEX-RK 3. Conversely, the implicit (EPIRK) solvers within IMEX-sEPIRK3(b) are up to 40 times faster at the same accuracy than the IMEX-RK 3 and the EPIRK scheme within IMEX-EPIRK is even up to 1200 times faster than the implicit Euler scheme within the IMEX-RK 1. Consequently, the key to efficiency indeed lies in the exponential integrators for the implicit treatment. Moreover, the experimental order of convergence for all IMEX schemes was verified for all three IMEX region decompositions (0%, 50%, 100% of implicit cells).
NACA0012 test case A steady state solution for the inviscid flow around the NACA0012 profile on a triangulation was successfully computed with a DG-P2 discretization and the new IMEX-sEPIRK3b scheme, see Fig. 3. The IMEX region was chosen to contain 50% of implicit cells. A CPU time comparison of the IMEX schemes until the time point of reaching a weighted residual of the density of 3.5 shows an efficiency gain of 54% of the IMEX-EPIRK scheme compared to IMEX-RK 1 and an efficiency gain of 28% of the IMEX-sEPIRK3b campared to the IMEX-RK 3 scheme.