1. Introduction
Highorder discretization methods enjoy increasing popularity because of their low dispersion errors and superior convergence properties [SE_Beck2014a, SE_Schaal2015a, TI_Tavelli2016a, SE_Winters2018a, SE_Fehn2019a, SE_Bassi2020a]. However, many studies use high order (3 to 16) in space, whereas time integration is often conducted with schemes of order 2, only. As the resolution of simulations increases, this approach requires extremely small time steps to match the temporal and spatial accuracy. Highorder timeintegration methods promise higher accuracy with larger steps and lower cost, which shall be assessed in the present paper.
While explicit methods offer the simplest approach, they become infeasible with growing spatial order due to the worsening stiffness of the diffusion terms [SE_Karniadakis2005a, SE_Canuto2011a]. This leaves the choice between implicit and semiimplicit methods, depending on the treatment of the convective terms. Stable implicit methods exist for virtually any order and possess the advantage of choosing the time step only according to accuracy or physical criteria [TI_Hairer1996a]. Several authors used Diagonally Implicit RungeKutta (DIRK) methods for simulation of compressible flows, e.g. [TI_Jothiprasad2003a, TI_Persson2010a, TI_Uranga2010a, TI_Pazner2017a, TI_Bassi2015a, TI_Pan2021a]. The application of implicit RungeKutta (IRK) methods to incompressible flows is more complicated due to the quasistatic nature of the continuity equation, but was investigated in [TI_John2006a, TI_Bassi2015a, TI_Noventa2016a]. In most of these studies, the order of convergence reached 3 to 5, although methods up to order 10 were considered in [TI_Southworth2021a]. As an alternative to IRK, variational methods such as Discontinuous Galerkin methods in time and spacetime element methods are gaining interest [TI_Tavelli2016a, TI_Ahmed2017a]. They use a polynomial expansion in time which makes them inherently implicit and capable to attain arbitrary order. A common drawback of implicit methods is the need to solve nonlinear systems in each step. These systems are linearized within a fixpoint or Newtontype method or by construction, as in the case of Rosenbrocktype RungeKutta methods [TI_John2006a, TI_Bassi2015a, TI_Noventa2016a]. Unfortunately, the resulting linear systems require computation and storage of the exact or approximate Jacobian of the nonlinear operator. Moreover, they are typically nonsymmetric and often illconditioned. Their fast solution depends on an efficient preconditioner, which may be hard to achieve. These difficulties diminish the benefits of implicit methods considerably.
Semiimplicit methods combine an implicit scheme for diffusion with an explicit one for convection and, hence, are often called IMEX methods. The explicit treatment of convection terms leads to a CourantFriedrichsLewy (CFL) condition, which depends on the type and the order of the spatial discretization. For methods based on polynomial elements, the admissible time step scales approximately as where is the element length, the fluid velocity and the polynomial degree, see e.g. [SE_Karniadakis2005a]. The implicit discretization of diffusive terms leads to symmetric linear systems, for which very efficient solution methods exist, e.g. [SE_Lottes2005a, SE_Janssen2011a, SE_Stiller2016b, SE_Stiller2017a]. This offers an enormous advantage over fully implicit methods, especially in high fidelity simulations, where the time step is limited by accuracy rather than stability constraints. Due to their simplicity, IMEX variants of linear multistep methods are popular, especially for incompressible flow simulation [TI_Karniadakis1991a, TI_Leriche2006a, TI_Klein2015a, TI_Fehn2017a]. Corresponding methods can be constructed for any order, however, increasing the latter beyond two leads results in the loss of Astability for the implicit part and worsens the stability restrictions for the explicit part [TI_Hairer1996a]. As an alternative, highorder IMEX RungeKutta (RK) methods have been proposed, e.g. [TI_Ascher1997a, TI_Kennedy2003a, TI_Nikitin2006a, TI_Cavaglieri2015a]. These methods have far better stability properties and a higher accuracy than the best IMEX multistep schemes of the same order. Several authors examined methods up to an order of 3 for compressible flow [TI_Kanevsky2007a, TI_Persson2011a], whereas applications of highorder RK methods to incompressible problems are rare. Unfortunately, the growing number of order conditions makes it difficult to construct IMEX RK methods of orders greater than 4. To the knowledge of the authors only one method of order 5 has been published [TI_Kennedy2003a]. In contrast, the spectral deferred correction (SDC) method can reach arbitrary convergence rates by incrementally increasing the order using sweeps of lower order schemes [TI_Dutt2000a, TI_Ong2020a]. TI_Minion2003b developed a semiimplicit spectral deferred correction (SISDC) method based on an IMEX Euler scheme. This method has excellent stability properties that are comparable with IMEX RK methods, but extends easily to higher order. It was adopted for incompressible flows in periodic domains in [TI_Minion2004a] and later extended to timedependent boundary conditions [TI_Minion2018a]. In [TI_Stiller2020a] the SISDC method was generalized to variable properties and demonstrated to reach convergence rates up to order 12 for threedimensional vortical flows with nonlinear (turbulent) viscosity.
The obvious shortcomings of loworder IMEX multistep methods lead to the question whether highorder RK and SDC methods can surpass and replace them as a suitable complement to higherorder spatial discretization methods. Facing this challenge, the present study investigates the convergence properties and computational efficiency of selected IMEX RK methods as standalone integrators and as predictors in SISDC methods for incompressible NavierStokes problems with constant or variable viscosity. To cope with the quasistatic incompressibility constraint and achieving a consistent treatment of nonlinear viscosity it is proposed to harnesses the explicit and implicit RK components as a partitioned RK method and to apply a tailored projection method on each stage. Spatial discretization is accomplished using a discontinuous Galerkin method of order 16. The convergence behavior of IMEX RK and SISDC methods up to order 6 is studied using various test cases featuring timedependent boundary conditions, nonlinear viscosity and turbulence. For each case an errorcost analysis is performed to assess the computational efficiency. A comparison with the widespread IMEX BDF2 method demonstrates the superiority of the proposed highorder timeintegration methods, even for direct simulations of turbulence with marginal spatial resolution.
The remainder of the paper is organized as follows: Section 2 describes the timeintegration methods using the example of a onedimensional convectiondiffusion equation with variable diffusivity. Section 3 extends these methods to incompressible NavierStokes problems, Section 4 presents the numerical experiments including a discussion of the results and Section 5 concludes the paper.
2. Timeintegration methods
2.1. Model problem
For explaining the IMEX approach and devising the semiimplicit treatment of nonlinear viscosity, the onedimensional convectiondiffusion equation is considered as a model problem:
(1) 
Here, the solution is a function in space and time, while the convection velocity and diffusivity may depend on . The simplest case with constant velocity and diffusivity possesses the spatially periodic solution with wave number . Its amplitude
satisfies the ordinary differential equation (ODE)
(2) 
where and represent the effect of convection and diffusion, respectively. The associated time scales and motivate a semiimplicit treatment as in the case of flow problems, since the decreases much faster with growing wave number than . Therefore, this problem serves as a model for characterizing the stability and accuracy of the investigated IMEX methods in Sec. 2.5.
A more general situation occurs, when a spatial discretization method is applied to the original problem (1
). In this case, the semidiscrete solution takes the form of a time dependent coefficient vector,
and is governed by the ODE system(3) 
where and are the discrete convection and diffusion operators, and and the coefficient vectors of velocity and diffusivity, both of which may depend on . As the time integration methods act individually on each component of (3), their description is based on the corresponding scalar problem
(4) 
where corresponds to and to , respectively. Note that the explicit dependence of on has been dropped for convenience.
2.2. IMEX methods for constant diffusivity
As a starting point, the timeintegration methods are introduced for constant diffusivity. This case allows for a clean IMEX approach, where the diffusive part is treated implicitly and the convective part explicitly. Using equidistant steps yields the times for which the approximate solution is sought. The IMEX backward difference formula of order 2 (IMEX BDF2) serves as a reference and is given by [TI_Frank1997a]
(5) 
where , , , and . Since the method is not selfstarting, it is initialized using an IMEX Euler step, i.e. and .
The IMEX RK methods investigated in this study are combinations of DIRK and explicit RungeKutta (ERK) schemes. They have an identical number of stages and different coefficients and , whereas the time nodes are synchronized, i.e. . These considerations yield the following general Butcher tableaux:
(6) 
To achieve a certain order, the DIRK and ERK methods must satisfy the corresponding order conditions as well as coupling conditions [TI_Kennedy2003a]. In addition to this, the IMEX RK methods selected for this study possess the firstsameaslast (FSAL) property, , and an Lstable DIRK part. The application to the model problem 4 with constant diffusivity is described in Algorithm 1. Table 1 provides an overview of the considered methods, including the number of stages, the theoretical order of convergence and further relations between the coefficients. The Butcher tableaux are given in Appendix A.
method  stages  order  

RKTR  3  2  
RKCB2  3  2  
RKCB3c  4  3  
RKCB3e  4  3  
RKCB4  6  4  
RKARS3  5  3 
RKTR is a 3stage version of the IMEX trapezoidal rule with better stability properties than the more common 2stage version. RKCB2, RKCB3c, RKCB3e and RKCB4 were developed by TI_Cavaglieri2015a, who named them IMEXRKCB2, IMEXRKCB3c, IMEXRKCB3e and IMEXRKCB4, respectively. The assumption of identical assembly coefficients reduced the number of coupling conditions and thus allowed to optimize these methods for stability and accuracy. Moreover, the DIRK part of RKCB4 attains stage order 2, which renders it less susceptible to order reduction when applied to stiff problems. RKARS3 was adopted from Ascher, Ruth and Spiteri [TI_Ascher1997a], who dubbed it (4,4,3), because the implicit and explicit part can be both reduced to 4 stages and the resulting IMEX method has order 3. As with RKTR, its final stage coincides with the assembly and hence eliminates the need for the latter. In applications to PDEs, like NavierStokes, this property allows for the exact fulfillment of boundary conditions, which is not guarantied for the assembled solution. Finally, it should be noted that all methods except RKCB3c satisfy and thus have a stiffly accurate DIRK part [TI_Hairer1996a].
2.3. Extension to variable diffusivity
The presence of a variable diffusivity renders the diffusion term nonlinear and, thus, complicates the solution of the implicit equations in each time step. With IMEX BDF2 this difficulty can be circumvented by extrapolating the diffusivity, i.e.
(7) 
where etc. This approach preserves second order accuracy and leads to discrete equations that are linear in and, hence, easier to solve as in the fully implicit case. As a possible side effect, diffusive stability restrictions are expected, but scarcely observed in practice.
Unfortunately, the semiimplicit treatment of nonlinear diffusion is not trivial with RungeKutta methods. Rosenbrocktype methods can handle variable coefficients, but are expensive due to the introduction of Jacobians and the more complicated structure of linear systems. TI_Boscarino2016a used partitioned RK methods to construct semiimplicit schemes without the need of nonlinear iterations. The IMEX RK method for variable diffusivity presented in Algorithm 2 was developed independently, but can be seen as a special case of this approach. At the begin of stage , it applies the ERK component to compute the preliminary solution (line 6), which is used to evaluate the diffusivity (line 7). The computation of the stage solution then resembles a fully implicit treatment of diffusion, but is actually semiimplicit, since is already known (line 8). As a consequence, each RK stage has a complexity comparable to the IMEX BDF2 method for variable diffusivity stated in (7), while the full scheme still retains the accuracy of the underlying IMEX RK method [TI_Boscarino2016a].
2.4. Spectral deferred correction
The SDC method divides the interval into subintervals with intermediate times . Following [TI_Dutt2000a] the GaussLobattoLegendre (GLL) points are chosen, such that and . A predictor sweep through all subintervals provides the initial approximations , which are joined in . The predictor can be any selfstarting timeintegration method [TI_Layton2007a]. In the present study, IMEX Euler and RK methods serve for this purpose. Starting from , the approximate solution is gradually improved by successive correction sweeps. Given the th approximation, , the corrector seeks a solution to the error equation [TI_Dutt2000a]
(8)  
where
is the Lagrange interpolant to
. Combining the equations for and leads to the correction equation(9)  
As proposed by TI_Minion2003b, the IMEX Euler rule is applied to the first integral, while the second one is approximated using with . Introducing and exploiting finally yields the discrete correction equation
(10)  
with subinterval length and weights , where is the Lagrange basis polynomial associated with . Each sweep solves Eq. (10) for to and increases the order by one [TI_Dutt2000a], although stiffness or timedependent boundary conditions can degrade convergence [TI_Boscarino2018a, TI_Minion2018a]. The converged solution satisfies the Lobatto IIIA collocation scheme with stages and, hence, achieves order at the end of the time interval [TI_Hagstrom2006a].
2.5. Linear stability and accuracy
For comparing the stability and accuracy of the considered timeintegration methods, the RHS of model problem (2) is split into an implicit and an explicit part such that
(11) 
The application of an IMEX onestep method yields the discrete evolution equation
(12) 
where is the (unknown) stability function, which depends on two complex parameters, and . Apart from the general case, the following special choices are possible:

implicit: , ,

explicit: , ,

semiimplicit: , .
Each of these cases yields a specific stability function that depends only on a single parameter, . Given a onestep IMEX method such as RK or SDC, the corresponding stability function can be evaluated for any by executing (12) once with , , and initial condition , such that . Additionally, this procedure provides the error . For a fair comparison of different methods the argument is scaled to the number of nontrivial substeps, i.e., for IMEX RK, for SISDC with Euler predictor and for SISDC with RK predictor. As is fixed for a given problem, the scaling applies also to the time step such that . Figure 3 shows the unscaled and scaled stability domains of selected methods for the semiimplicit case and . The threshold is of interest because the imaginary value
corresponds to the eigenvaluebased CourantFriedrichsLewy number, i.e.,
. Its critical value depends on the mangnitude of real part , which is sometimes referred to as Diffusion number.As expected, the sixthorder SDC method admits the larges time steps over a wide range of (Fig. (a)a). However, when applying the substep scaling (Fig. (b)b), RKCB3e allows scaled time steps roughly 1.5 times larger than RKARS3 and 4 times larger than SDCEu(3,5). SDCCB3e(3,3) and SDCARS3(3,3) resemble the latter and are omitted for clarity. All remaining RK methods, including CB2, range between ARS3 and CB3e. Figure 6 compares the scaled accuracy domains of the same methods for two different thresholds. For errors smaller than , RKCB3e is clearly the most efficient method, followed by RKCB4 and then RKCB3c and SDCEu(3,5), with considerably smaller or, respectively, . With increasing accuracy, the 6th order SDC methods gain efficiency and become competitive with RKCB3e for (Fig. (b)b). In summary, its superior stability and accuracy render RKCB3e the most promising method for the considered model problem. However, this clear dominance is not necessarily preserved in the application to incompressible flow problems, where nonlinearity, boundary conditions, splitting errors and spatial discretization introduce additional challenges that can affect stability as well as accuracy of the methods.
3. Application to incompressible NavierStokes problems
3.1. Governing equations
The timeintegration methods introduced above are applied to incompressible flows with constant density and variable viscosity in a simply connected spatial domain . The velocity satisfies the momentum (NavierStokes) and continuity equations
(13)  
(14) 
in , where is the pressure and
the viscous stress tensor, both divided by density;
the kinematic viscosity, which may depend on , and a known forcing term. The flow problem is closed by initial and boundary conditions(15)  
(16) 
satisfying the continuity and compatibility constraints, and , respectively.
3.2. Time integration
The IMEX methods are adapted to incompressible NavierStokes problems by applying a projectionbased splitting scheme in each step and with RK at each stage. This scheme is based on the rotational velocity correction scheme developed in [TI_Guermond2003b] and extended to variable viscosity in [TI_Stiller2020a]. For stating the methods, the flow equations are rewritten in the form
(17) 
where with , and
(18a)  
(18b) 
The additional viscous contribution introduces a divergence penalty which, for constant viscosity, yields the laplacian form of the viscous term, i.e., . Semidiscrete quantities write as follows in general, with RK at stage and with SDC at time after corrections. Intermediate results are indicated by one or more primes, e.g. .
3.2.1. Imex Bdf2
With IMEX BDF2, the splitting scheme is applied once per time step, which yields four substeps:

Extrapolation
(19) 
Projection
(20) 
Diffusion
(21) 
Additional projection
(22)
The extrapolation (19) can be viewed as an explicit step that includes all terms except forcing, which is treated implicitly. Adding yields an extra divergence penalty and transforms the diffusion term into rotational form, , if is constant. The projection step (20) requires the pressure , which is obtained by solving the related Poisson problem
(23) 
Application of in the first part of (20) yields the divergencefree intermediate velocity . The following step solves the implicit diffusion problem (21), which removes the explicit approximation of and the extra divergence penalty introduced in (19). If is constant, the diffusion step preserves continuity. However, a variable viscosity can result in a nonsolenoidal contribution to the RHS in (21) and thus affect the continuity of the velocity field . Therefore, it is generally necessary to conclude the time step with an additional projection according to Eq. (22). If the viscosity is constant, the projection is skipped such that .
3.2.2. Imex Rk
Algorithm 3 outlines the IMEX RungeKutta method using a splitting similar to BDF2 at each stage except the first. In comparison to Algorithm 2, takes the role of and is used to compute the viscosity . The auxiliary quantity combines the pressure contributions of different stages. This simplifies the handling and removes the need for keeping the stagelevel pressures. Since no initial value is required, the computation of is dropped as well. Line 11 of the algorithm presents an incremental form of the assembly which allows to omit terms with vanishing factors . The final projection (line 13–14) removes the residual velocity divergence in case of variable viscosity and is skipped with constant .
3.2.3. Sisdc
The spectral deferred correction method considered in this study uses either IMEX Euler or a RungeKutta method as the predictor. Each time step is divided into subintervals such that the intermediate times represent the GLL points, as described in Sec. 2.4. Starting from , the predictor sweeps through the subintervals and computes for to . This initial approximation is improved by the IMEX Euler corrector outlined in Algorithm 4. It is a streamlined version of the most promising approach devised in [TI_Stiller2020a] with simplified pressure handling analogous to IMEX RK. Except for the extrapolation that contains additional SDC terms (line 17), a corrector step over one subinterval corresponds to an ordinary IMEX Euler step, performing the projection on lines 1819, followed by the diffusion step on line 20 and, in case of variable , an additional projection on lines 2223. After sweeps, the corrector returns the solution and thus completes the time step.
Comments
There are no comments yet.