Assessment of high-order IMEX methods for incompressible flow

by   Montadhar Guesmi, et al.
TU Dresden

This paper investigates the competitiveness of semi-implicit Runge-Kutta (RK) and spectral deferred correction (SDC) time-integration methods up to order six for incompressible Navier-Stokes problems in conjunction with a high-order discontinuous Galerkin method for space discretization. It is proposed to harness the implicit and explicit RK parts as a partitioned scheme, which provides a natural basis for the underlying projection scheme and yields a straight-forward approach for accommodating nonlinear viscosity. Numerical experiments on laminar flow, variable viscosity and transition to turbulence are carried out to assess accuracy, convergence and computational efficiency. Although the methods of order 3 or higher are susceptible to order reduction due to time-dependent boundary conditions, two third-order RK methods are identified that perform well in all test cases and clearly surpass all second-order schemes including the popular extrapolated backward difference method. The considered SDC methods are more accurate than the RK methods, but become competitive only for relative errors smaller than ca 10^-5.



There are no comments yet.


page 16


A spectral deferred correction method for incompressible flow with variable viscosity

This paper presents a semi-implicit spectral deferred correction (SDC) m...

High-Order Discretization of Backward Anisotropic Diffusion and Application to Image Processing

Anisotropic diffusion is a well recognized tool in digital image process...

Fast multigrid solution of high-order accurate multi-phase Stokes problems

A fast multigrid solver is presented for high-order accurate Stokes prob...

Implicit Methods with Reduced Memory for Thermal Radiative Transfer

This paper presents approximation methods for time-dependent thermal rad...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

High-order 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. High-order time-integration 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 semi-implicit 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 Runge-Kutta (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 Runge-Kutta (IRK) methods to incompressible flows is more complicated due to the quasi-static 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 space-time 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 fix-point or Newton-type method or by construction, as in the case of Rosenbrock-type Runge-Kutta 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 non-symmetric and often ill-conditioned. Their fast solution depends on an efficient preconditioner, which may be hard to achieve. These difficulties diminish the benefits of implicit methods considerably.

Semi-implicit 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 Courant-Friedrichs-Lewy (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 A-stability for the implicit part and worsens the stability restrictions for the explicit part [TI_Hairer1996a]. As an alternative, high-order IMEX Runge-Kutta (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 high-order 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 semi-implicit 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 time-dependent 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 three-dimensional vortical flows with nonlinear (turbulent) viscosity.

The obvious shortcomings of low-order IMEX multistep methods lead to the question whether high-order RK and SDC methods can surpass and replace them as a suitable complement to higher-order spatial discretization methods. Facing this challenge, the present study investigates the convergence properties and computational efficiency of selected IMEX RK methods as stand-alone integrators and as predictors in SISDC methods for incompressible Navier-Stokes problems with constant or variable viscosity. To cope with the quasi-static 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 time-dependent boundary conditions, nonlinear viscosity and turbulence. For each case an error-cost analysis is performed to assess the computational efficiency. A comparison with the widespread IMEX BDF2 method demonstrates the superiority of the proposed high-order time-integration methods, even for direct simulations of turbulence with marginal spatial resolution.

The remainder of the paper is organized as follows: Section 2 describes the time-integration methods using the example of a one-dimensional convection-diffusion equation with variable diffusivity. Section 3 extends these methods to incompressible Navier-Stokes problems, Section 4 presents the numerical experiments including a discussion of the results and Section 5 concludes the paper.

2. Time-integration methods

2.1. Model problem

For explaining the IMEX approach and devising the semi-implicit treatment of nonlinear viscosity, the one-dimensional convection-diffusion equation is considered as a model problem:


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)


where and represent the effect of convection and diffusion, respectively. The associated time scales and motivate a semi-implicit 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 semi-discrete solution takes the form of a time dependent coefficient vector,

and is governed by the ODE system


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


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 time-integration 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]


where , , , and . Since the method is not self-starting, 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 Runge-Kutta (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:


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 first-same-as-last (FSAL) property, , and an L-stable 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
RK-TR 3 2
RK-CB2 3 2
RK-CB3c 4 3
RK-CB3e 4 3
RK-CB4 6 4
RK-ARS3 5 3
Table 1. Properties of IMEX Runge-Kutta methods.

RK-TR is a 3-stage version of the IMEX trapezoidal rule with better stability properties than the more common 2-stage version. RK-CB2, RK-CB3c, RK-CB3e and RK-CB4 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 RK-CB4 attains stage order 2, which renders it less susceptible to order reduction when applied to stiff problems. RK-ARS3 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 RK-TR, its final stage coincides with the assembly and hence eliminates the need for the latter. In applications to PDEs, like Navier-Stokes, 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 RK-CB3c satisfy and thus have a stiffly accurate DIRK part [TI_Hairer1996a].

1:procedure ImexRK()
3:     for  do
5:     end for
7:end procedure
Algorithm 1 IMEX Runge-Kutta method for model problem with constant diffusivity.

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.


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 semi-implicit treatment of nonlinear diffusion is not trivial with Runge-Kutta methods. Rosenbrock-type 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 semi-implicit 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 semi-implicit, 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].

1:procedure ImexRK()
5:     for  do
9:     end for
11:end procedure
Algorithm 2 Semi-implicit partitioned Runge-Kutta method for variable diffusivity.

2.4. Spectral deferred correction

The SDC method divides the interval into subintervals with intermediate times . Following [TI_Dutt2000a] the Gauss-Lobatto-Legendre (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 self-starting time-integration 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]



is the Lagrange interpolant to

. Combining the equations for and leads to the correction equation


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


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 time-dependent 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 time-integration methods, the RHS of model problem (2) is split into an implicit and an explicit part such that


The application of an IMEX one-step method yields the discrete evolution equation


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: , ,

  • semi-implicit: , .

Each of these cases yields a specific stability function that depends only on a single parameter, . Given a one-step 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 semi-implicit case and . The threshold is of interest because the imaginary value

corresponds to the eigenvalue-based Courant-Friedrichs-Lewy number, i.e.,

. Its critical value depends on the mangnitude of real part , which is sometimes referred to as Diffusion number.

As expected, the sixth-order SDC method admits the larges time steps over a wide range of (Fig. (a)a). However, when applying the substep scaling (Fig. (b)b), RK-CB3e allows scaled time steps roughly 1.5 times larger than RK-ARS3 and 4 times larger than SDC-Eu(3,5). SDC-CB3e(3,3) and SDC-ARS3(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 , RK-CB3e is clearly the most efficient method, followed by RK-CB4 and then RK-CB3c and SDC-Eu(3,5), with considerably smaller or, respectively, . With increasing accuracy, the 6-th order SDC methods gain efficiency and become competitive with RK-CB3e for (Fig. (b)b). In summary, its superior stability and accuracy render RK-CB3e 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.

Figure 3. Unscaled and scaled stability domains of selected semi-implicit methods.
Figure 6. Scaled accuracy domains of selected semi-implicit methods.

3. Application to incompressible Navier-Stokes problems

3.1. Governing equations

The time-integration 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 (Navier-Stokes) and continuity equations


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


satisfying the continuity and compatibility constraints, and , respectively.

3.2. Time integration

The IMEX methods are adapted to incompressible Navier-Stokes problems by applying a projection-based 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


where with , and


The additional viscous contribution introduces a divergence penalty which, for constant viscosity, yields the laplacian form of the viscous term, i.e., . Semi-discrete 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:

  1. Extrapolation

  2. Projection

  3. Diffusion

  4. Additional projection


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


Application of in the first part of (20) yields the divergence-free 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 non-solenoidal 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 Runge-Kutta 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 stage-level 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 .

1:procedure ImexRK()
3:     for  do
10:     end for
12:     if final projection then
15:     end if
16:end procedure
Algorithm 3 IMEX Runge-Kutta method for incompressible flow with variable viscosity. Compact notation is applied were possible, e.g. .

3.2.3. Sisdc

The spectral deferred correction method considered in this study uses either IMEX Euler or a Runge-Kutta 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 18-19, followed by the diffusion step on line 20 and, in case of variable , an additional projection on lines 22-23. After sweeps, the corrector returns the solution and thus completes the time step.

1:procedure SISDC()
2:     Initialize and
4:     for  do initial viscosity and RHS
6:         if  then 
7:         if  then 
8:     end for
9:     for  do correction sweeps
13:         for  do sweep through subintervals
21:              if final projection then