Development of a Balanced Adaptive Time-Stepping Strategy Based on an Implicit JFNK-DG Compressible Flow Solver

A balanced adaptive time-stepping strategy is implemented in an implicit discontinuous Galerkin solver to guarantee the temporal accuracy of unsteady simulations. A proper relation between the spatial, temporal and iterative errors generated within one time step is constructed. With an estimate of temporal and spatial error using an embedded Runge-Kutta scheme and a higher order spatial discretization, an adaptive time-stepping strategy is proposed based on the idea that the time step should be the maximum without obviously influencing the total error of the discretization. The designed adaptive time-stepping strategy is then tested in various types of problems including isentropic vortex convection, steady-state flow past a flat plate, Taylor-Green vortex and turbulent flow over a circular cylinder at Re=3900\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Re}=3\,900$$\end{document}. The results indicate that the adaptive time-stepping strategy can maintain that the discretization error is dominated by the spatial error and relatively high efficiency is obtained for unsteady and steady, well-resolved and under-resolved simulations.


Introduction
Implicit time integration schemes have the potential to increase the efficiency of high-Reynolds number simulations by relaxing the challenging time step restriction, and have been successfully adopted in a variety of steady-state and unsteady simulations [2,21,22,27,29]. However, implicit time integration schemes are generally more complex and introduce many parameters, which usually have large influences on the performance of the solver. Take the widely used implicit time integration schemes based on Runge-Kutta (RK) temporal discretization scheme and Jacobian-free Newton Krylov (JFNK) method as an example, there are parameters such as the order of RK scheme, the time step size, the Newton iteration convergence tolerance, the linear iteration convergence tolerance and possibly other parameters introduced by preconditioners [15,16,29]. These parameters are usually highly problem dependent and have large influences on the accuracy, efficiency and robustness in specific simulations, the choices of which actually become a complex multi-objective optimization problem. Therefore, a reliable approach for determining these parameters is essential especially for computational fluid dynamics (CFD) software with a wide range of application areas and for the automatic simulations of massive cases in a design pipeline.
In steady-state simulations, there were many studies on this topic, most of which mainly focus on accelerating the convergence to steady state. In reference [28], an expert system for choosing an efficient CFL was developed, in which the convergence history is manually separated into three different stages and different adaptive strategy is used in each stage. In reference [18], a solution-limiting method is developed to determine the CFL. Bücker et al. compared some of the existing adaptive methods in 2009 and concluded that there is no clear winner [6]. In 2019, an error based adaptive CFL method was developed to accelerate the late stage of convergence but still need to manually divide the convergence history into different stages [12]. Compared with the studies on efficiency, much fewer studies have focused on robustness. Lian et al. [18] addressed this problem with the solution-limiting method and claimed that the robustness is largely improved. Other attempts on this topic are mainly some techniques such as rolling back and recomputing with a smaller time step if the solution is not satisfactory [30]. These studies have highlighted the importance of these parameter choices on different aspects of the solver and on different research areas. Although, significant progress has been made, as is concluded in reference [6] optimal CFL evolution is still an open problem.
For unsteady simulations, there are more parameters and there is an extra trouble to worry about, the temporal accuracy. The idea of using adaptive time-stepping in unsteady simulations is old, which can date back to time adaptation studies of ordinary differential equations (ODEs) [1,24,25]. Various methods for choosing the step size (or step size controllers) were proposed and compared in the history as well as other topics such as the requirement of the temporal error estimator and the stability properties in stiff problems, a comprehensive review of which can be found in reference [24,25]. Birken adopted a similar idea as those in the ODE systems and developed a time step adaptation method in the simulation of compressible Navier-Stokes equations using an embedded scheme to estimate the local temporal error. This method is adopted in the comparison of Rosenbrock methods and explicit first stage singly diagonally implicit Runge-Kutta (ESDIRK) time integration methods [5]. Noventa et al. developed a similar time step adaptation strategy in unsteady incompressible turbulent flow simulations [21], which is further tested in compressible simulations focusing on the comparison of different step size controllers [22]. The adaptive method based on local temporal error was further developed for goal oriented time step adaptation and compared with a method based on global error estimate [19]. However, most of them such as, the studies in references [4,5,21,22], consider the temporal discretization methods separately as an ODE solver using the method of lines but do not take into account specific physical properties of the problem and specific spatial discretization adopted. These methods rely on a user defined temporal error tolerance, the choice of which is highly problem dependent and needs a deep understanding of the simulation settings. However, a naive choice of temporal error tolerance may lead to a large increase of CPU time without obviously improving the results as has been demonstrated in reference [21].
In this work, an adaptive time-stepping method is proposed to maintain temporal accuracy and high efficiency. The basic idea is to determine the time step based on the relations of different discretization errors. The time step is determined as the largest time step such that temporal errors are smaller than the spatial discretization errors. In this sense, the temporal discretization is believed to be sufficiently accurate since further decreasing the time step will not improve the simulation obviously. As will be shown in numerical tests in Sect. 3, this choice is the relatively efficient choice because it is around the maximum choice to be sufficiently accurate in time in our definition.
The governing equations, numerical methods and numerical error relations are presented in Sect. 2. Moreover, the adaptive strategy proposed in this paper is also presented in Sect. 2. Results of various test cases are presented in Sect. 3 to demonstrate the properties of the adaptive methods in different problems. Finally, conclusions are drawn in Sect. 4.

Governing Equations and Numerical Methods
The governing compressible Navier-Stokes equations, representing conservation of mass, momentum and energy, can be written in abridged form as where L is the analytical nonlinear spatial operator, d is the number of dimensions of the problem, and = [ , u 1 , ⋯ , u d , E] T is the vector of conservative variables. The inviscid flux, = ( ) , and the viscous flux, = ( , ∇ ) , in the ith direction are given by

3
Here, E = (C v T + u i u i ∕2) is the total energy in which C v is the specific heat capability with constant volume, p = ( − 1)(E − u i u i ∕2) is the pressure, is the specific heat ratio, is the heat flux, is the dynamic viscosity, and is the thermal conductivity.
In what follows we briefly outline the formulation of a compressible flow solver based on a discontinuous Galerkin (DG) spatial discretization [13], an implicit time integration via ESDIRK temporal discretization schemes [17] and a JFNK method [16]. Further details of the solver can be found in reference [29].

Spatial Discretization
In DG methods, the computational domain is divided into non-overlapping elements and a polynomial space is defined on each element. Performing integration by parts after multiplying the test function p to Eq. (1), we get the weak form of Eq. (1) in element e: where e is the element boundary, and n is the boundary flux on the elemental outward normal direction ( ). Furthermore, the vector of conserved variables is expressed in the polynomial space as on each computational element. Here, p is the coefficient of the pth basis function and N P + 1 is the number of basis functions. To discretize Eq. (3), a quadrature rule is adopted to approximate the integration. The on these quadrature points can be calculated by where the subscript i indicates values on the ith quadrature point, and is the backward transform matrix. Using the Roe's Riemann flux [26] and the symmetric interior penalty DG (SIPG) method [10] to approximate the advective numerical flux and the viscous numerical fluxes, respectively, the semi-discrete equation can be expressed as where w is the quadrature weight, J is the grid metric Jacobian, represents a diagonal matrix, = T (wJ) is the mass matrix, j is the derivative matrix in the jth direction, the superscript represents the corresponding matrices or variables on element boundaries, ̂ n j is a symmetric flux from the SIPG method, is the mapping matrix between and , and is the interpolation matrix from quadrature points of an element to quadrature points of its element boundaries. Details of the spatial discretization can be found in reference [29].
The corresponding semi-discrete equation for the vector can be obtained using relation of Eq. (5) and which is written as The spatial truncation error can be expressed as where C s is a coefficient related to the discretization scheme but independent of , D s ( ) is a spatial derivative term dependent on and P is the polynomial order of the basis functions in Eq. (5).

Temporal Discretization
An embedded version of the ESDIRK method [17] is adopted to discretize Eq. (3) in time. The ith stages of the process of ESDIRK are where a ij , b i and b i are coefficients of the Butcher tableaus on the ith RK stage which can be found in Appendix A and in reference [17]. n+1 is the approximation at time step n + 1 . ̂ n+1 is the embedded solution for temporal error estimate at the cost of one extra implicit stage. For the embedded ESDIRK schemes adopted, ̂ n+1 is one order higher than n+1 . Thus, the leading term of the temporal error of n+1 can be accurately estimated as where b S+1 is always zero in the ESDIRK schemes adopted, C t is a coefficient related to the ESDIRK scheme adopted but independent of , D t ( ) is a problem dependent temporal derivative term, and N is the order of accuracy of the ESDIRK scheme. In the derivation of Eq. (15), the analytical equations, Eq. (1), have been adopted to replace the spatial operator with the temporal derivatives, thus spatial error is not included. The estimation also excludes the temporal errors from previous time steps and will be referred to as the local temporal error in this paper.

Discrete Formulation
The discrete formulation of the governing equations is obtained by replacing and L in Eqs. (10)- (14) with their discrete counterparts ( and L ). To reduce the size of the nonlinear system for implicit stages ( a ii ≠ 0 ), we choose to rewrite the discrete counterparts of Eqs. (11)-(12) as since the length of is usually smaller than that of . In rewriting the above equations, the relations of Eqs. (5) and (7) are used.

JFNK Method
Equation (17) of the implicit stages ( a ii ≠ 0 ) can be written as the nonlinear system This system is iteratively solved by Newton's method with an initial guess 0 = (i) . The Newton step, Δ = k+1 − k , is updated through the solution of the linear system A preconditioned GMRES iterative solver is adopted here to solve the unsymmetric system of Eq. (20). A Jacobian-free method is further adopted in the GMRES to approximate the Jacobian matrix and a vector inner product operator using a finite difference approximation, that is where i is a set of orthogonal vectors in the Krylov space, the expression of can be found in references [16,29]. A low memory block relaxed Jacobi iterative preconditioner is used in the GMRES solver to speedup the simulation [29]. The GMRES tolerance of 0.05‖ ( k )‖ adopted here has no influence on the temporal accuracy but affects its efficiency. Other choices of the GMRES tolerance can be found in references [8,14].
When the Newton residual is smaller than the convergence tolerance of the Newton iterations , which is (i) it = k is regarded as the approximate solution of the nonlinear system, where is the remaining Newton residual vector ( (i) it ) after convergence. A proper choice of is important because it determines the magnitude of the iterative error introduced by the iterative solver, which may degrade the overall error and the convergence order of accuracy [22]. Assume (i) is the exact solution of the nonlinear system and (i) it is the iterative approximation, we can get ( (i) ) = and ( (i) it ) = (i) . Thus, The contribution of the iterative error to the total discretization error will be discussed in the following.

Error Estimate
In this section, the local discretization errors introduced at time step n + 1 are estimated, which mainly consist of the spatial error, temporal error and iterative error. The accumulation of these errors, which we refer to as global errors, is more difficult to estimate and will not be considered in this paper. Studies considering the global errors can be found in reference [19]. We recall that the numerical solution at the time step n + 1 is calculated by (20) and the exact solution at n + 1 , (t n+1 ) , based on n can be expressed as where O(Δt N+2 ) represents all the high-order truncation terms and is not detailed here. Subtracting Eq. (24) by Eq. (25) and ignoring the high-order truncation terms of O(Δt N+2 ) , the leading terms of the local discretization error can be expressed as If we now add and subtract the summation terms of Δt to the right-hand side of Eq. (26), we obtain where b S+1 is always zero in the ESDIRK schemes adopted. The L  (28) suggests that the error generated in one time step is the sum of the local temporal error, averaged spatial error and averaged iterative error. This indicates that to decrease the total error one should decrease these three errors simultaneously and it will not help much if one of them is decreased much smaller than the largest of the three terms.
Based on this observation, an adaptive time-stepping strategy and an adaptive Newton tolerance are developed as follows in Sects. 2.2 and 2.3, respectively. The basic idea is to choose a time step and a Newton tolerance such that the temporal and iterative errors are smaller than the spatial error. In Sect. 2.2, the iterative error is assumed to be negligible compared with the spatial and temporal errors, which can be achieved as long as Eq. (19) converges to a sufficiently small level. The method to control the iterative error will be discussed in Sect. 2.3.

Adaptive Time-Stepping Strategy
Based on the error estimate, a time step adaptation method is proposed with the basic idea that the time step should be the maximum value without obviously influencing the total discretization error. With the spatial error fixed in a specific simulation, this implies that with < 1.0 and a proper defined norm. When ‖ n+1 t ‖ is orders of magnitude smaller than Δt‖̄ s ‖ , the total error will be dominated by the spatial error and further decreasing the time step (thus the temporal error) will not improve the solution much.
To maintain Eq. (29) is satisfied, spatial and temporal errors need to be estimated. The embedded ESDIRK scheme offers an accurate estimate of temporal error n+1 t using Eq. (15). As for the spatial truncation error ̄ s , it is estimated using a one order higher approximation at the end of each time step, that is the leading term of which is the same as that in Eq. (9). Equation (29) implies that the desired time step ( Δt ) should satisfy that Dividing Eq. (31) by Eq. (15), time step Δt can be calculated by which can be used as the time step for the next time step. The simplest elementary controller [25] is adopted here in determining the time step and other choices are also possible [25]. In the derivation of Eqs. (30), (31) and (32), the change of D t ( ) is assumed to be small throughout the step.
In the designs of adaptive time-stepping methods in references [5,22], an L 2 norm of the n+1 t for all the variables throughout the flow field is adopted as the norm in their adaptive strategy. However, the norm of the error with this choice may be dominated by a specific variable or the large scale unsteady structures. As a result, the adapted time step may be determined by one of the variables and/or the dominating unsteady structure in the flow field, which could be too large for the small unsteady structures of interest.
To avoid these problems, the norm is defined in each element. The L 2 norms of n+1 t and ̄ s are calculated in element e for variable m and the elemental time step is calculated using Here, the m is added to the denominator to avoid division by zero. But it also serves as a threshold of the error magnitude of interest. If the local temporal error magnitudes are much smaller than m , the local unsteady structures are either too weak to be of interest, or are already accurately captured, which indicates the current time step choice is already . excessively accurate in that element. Thus, a 1.5 N m is added to the numerator to maintain the elemental and the time step will at least grow by 1.5 times. The choice of 1.5, which is not optimized, ensures that the time step will increase but not too fast.
The value of m adopted here is where rms it,m denotes the root mean square value over the whole flow field of the mass, momentum magnitude and energy variables. Such a threshold value is based on the rounding error of the mth variable. One purpose of introducing this term is to avoid noisy elemental time steps caused by rounding errors in uniform flow areas.
With elemental time steps Δt n+1 e,m , the time step for variable m can be calculated using different strategies. For example, a minimum of Δt n+1 e,m would maintain that Eq. (29) is satisfied in every element. Here we choose to calculate it by a weighted average of Δt n+1 e,m as The temporal error ‖ t ‖ r e,m is used as the weights such that elements with larger temporal errors have larger contributions to the global time step. The r ∈ [0, ∞) can be chosen based on the smallest temporal structure of interest in the simulations, a larger value of which would lead to a larger contribution from elements with large temporal errors. In this paper, r = 1 is used. This weighted averaging helps maintain a relatively smooth change of the time step. At last, the time step Δt n+1 is calculated by

Adaptive Newton Tolerance
A similar idea as the time step adaptation can be adopted here in choosing the Newton tolerance so that the iterative error is smaller than temporal error to avoid affecting the temporal accuracy. In references [5,21], a Newton tolerance in Eq. (22) is proposed by directly relating the Newton tolerance to the temporal error, that has the form with = 0.1 . However, the relation between the Newton residual vector, , and the iteration error, Δt̄ it , introduced into Eq. (28) is not trivial. In the following this relation is estimated to illustrate the assumptions made in choosing the Newton tolerance of Eq. (37).
which maintains that the iterative error is smaller than the temporal error in Eq. (28)  The presence of L � ∕Δt − a ii L � −1 , which is difficult or expensive to estimate accurately, makes the relation between Δt it and problem dependent and difficult to evaluate. As Δt becomes larger L � ∕Δt − a ii L � −1 becomes closer and closer to a diagonal matrix with diagonal values of 1∕a ii and at the limit of Δt → ∞ , Eq. (44) becomes with 0.25 < a ii < 0.45 in the ESDIRK schemes adopted. When Δt becomes smaller and smaller, the entities of the matrix L � ∕Δt − a ii L � −1 will tend to zeros. Therefore, an even larger Newton tolerance will maintain that Eq. (44) is satisfied. Overall, Eq. (45) appears to provide a good estimate of the convergence upper bound at least at the limits of Δt → ∞ and Δt → 0 . Based on numerical experiments we observe that the Newton tolerance of Eq. (37) seems good enough to maintain temporal convergence accuracy, as will be shown in Sect. 3.1. However, the assumptions made in the derivation of the Newton tolerance may be invalid in some simulations, and if accurate estimate of Newton tolerance are needed, L � ( ∕Δt − a ii L � ) −1 can be estimated using the Jacobian matrix or simply by a calibration for a specific problem.
Although based on some assumptions, Eq. (37) is still much less problem dependent compared with Newton tolerances commonly used such as (43) as is shown in references [5,21,22] and will be shown in Sect. 3.1. It will also change in a more consistent manner with the temporal error when Δt changes. For simulations not using the adaptive Newton tolerance in the following, 2 = 10 −3 in Eq. (46) is adopted unless otherwise specified.

Summary of the Adaptive Strategy and Discussion
The basic idea behind the time step adaptation strategy proposed in Sect. 2.2 is the relation between temporal and spatial errors that is sufficient for maintaining temporal accuracy in unsteady simulations. Here, we give a brief discussion to illustrate that it is a rational choice because it shares a very similar logic behind the convergence tests in CFD, which is also a key point in maintaining the high efficiency of implicit time integration methods.
The implicit time integration methods could be more efficient that the explicit methods in some stiff problems such as problems with low Mach number, high Reynolds number and highly stretched grids [2,21,22,27,29]. The assumption behind this conclusion is that we can safely increase the time steps and temporal errors when using implicit time integration methods without obviously degrading the simulation results. In another word, the comparison of implicit and explicit time integration in efficiency is not based on the computation time of the same temporal error level as is done in most studies of time integrations methods [3,11] but is based on the computation time as long as the results are sufficiently accurate 1 . This implicitly defined sufficiently accurate condition is the maximum error level without obviously degrading the target quantity of interest, which is the one we are seeking when performing convergence tests in CFD simulations. However, this condition usually depends on the unsteady properties of the problem (the magnitude and frequency of the unsteady waves), the scale of interest, the acceptable error level, etc, and is not known before an accurate simulation result is obtained.
Similarly, the error condition behind the time adaptation strategy in Sect. 2.2, is the maximum temporal error level without obviously degrading the results of the discrete partial differential equation (PDE) system. Although, convergence tests are still needed to maintain the accuracy physically, the strategy avoids improper error relations, which may seriously lower the efficiency without obviously improving the results. In a case presented in reference [21], the improper relations of the spatial and temporal errors can increase the CPU time by about 20 times without obviously improving the results.
In some adaptive time-stepping strategies, the adapted time step is limited to be no larger than a time step Δt limit , which is user defined based on their understanding (46) of the physical time scale to be resolved as in reference [22]. However, Eqs. (28) and (29) indicate that limiting the time step will not improve the results but possibly will increase the computational costs. Instead, refining the spatial grid will be much more efficient in improving the results. Therefore, it is advised to refine the mesh instead of limiting the time step to improve the computation accuracy in our time adapting strategy when the adaptive time step is obviously too large to capture the unsteady flow.
As a by-product of this observation, the adapted time step can also serve as a good indicator of mesh refinement if the temporal scale of interest is well defined and the properties of the ODE solver (ESDIRK in this paper) are well understood. For a specific ODE solver, such as DIRK, we can assess the ODE solver in predicting a single wave at a specific time with a given time step [11]. If the unsteady flow scale (unsteady frequency) of interest is clearly defined in our simulations, this allows us to determine if the unsteady flow scale is satisfactorily captured with the ODE solver and the adapted Δt . If the temporal error is too large to be satisfactory, the relation (29) tells us that decreasing the time step cannot improve the results, and that a more efficient way to increase accuracy is to refine the spatial mesh instead. Therefore, we can use the adaptive time steps as a criteria for spatial mesh refinement.
In the adaptative strategy, extra computational costs are required to estimate spatial and temporal errors. The embedded scheme of Eq. (14) requires an extra implicit RK stage, which will approximately increase the computational cost by 1∕(S − 1) , where S is the number of stages of the main ESDIRK scheme with values of 3, 4 and 6 for the second, third and fourth order ESDIRK schemes adopted, respectively. Similarly, the spatial error estimation in Eq. (30) requires the calculation of a higher order spatial discretization operator L P+1 n+1 it . The extra cost from spatial error estimation is relatively lower than those of temporal error estimation for stiff problems since up to an order of magnitude or more L operations are needed in one implicit RK stage. For example, in simulations of P = 2 and P = 3 , the extra cost of L P+1 n+1 it is only about 1.5 times the cost of L P n+1 it , while the extra cost in temporal error estimating is more 10 times of that. Therefore, we can roughly estimate the extra cost to be 100∕(S − 1)% . For simulations with little or slow time scale changes, a freezing of the adapted time step can be adopted to lower the additional costs.
Compared with adaptive time-stepping strategies in references [4,5,21,22], the current strategy is different in the following ways.
• The adaptive time step is determined by the relation between the temporal and spatial errors instead of a user defined temporal error level. Therefore, highly problem dependent user defined tolerances are avoided. • An embedded ESDIRK scheme with an order of N + 1 is employed instead of schemes of (N − 1) th order, which predicts the temporal error more accurately. This also helps avoid problem dependent calibration parameters of the estimated errors. • The adaptation is determined element by element first instead of using the L 2 norm of the temporal error of the whole flow field. Therefore, it has the potential to avoid time step being determined primarily by the large scale unsteady structures and avoid smallscale unsteady structures being under-resolved.

Numerical Results
Numerical experiments are conducted in this section to demonstrate the properties of the methods proposed in Sect. 2. The isentropic vortex is simulated in Sect. 3.1 to validate the Newton tolerance choice and illustrate some basic properties of the adaptive time-stepping. The steady-state flow past a flat plate is adopted in Sect. 3.2 to demonstrate the performance of the adaptive method in steady-state simulations. In Sects. 3.3 and 3.4, the Taylor-Green vortex problem and the turbulent flow over a circular cylinder at Re = 3 900 are studied to illustrate the performance of the methods in freely decaying and wall-bounded turbulent flow simulations, respectively.

Isentropic Vortex Convection
This test case is adopted because its analytical solution is available, which makes it easier for spatial and temporal error estimations. In the computational domain 16 π 2 e 2(1−r 2 ) to the boundaries in both directions. The uniform mesh and the density distribution are illustrated in Fig. 1.
To estimate different errors, different solutions are adopted as the reference in the error estimation. The total discretization error of a solution is obtained by subtracting it with the analytical solution of Eq. (47). As in reference [3], the difference between the solution and a solution with a very small time step is regarded as the sum of temporal and iterative errors as the two solutions share the same spatial discretization. The difference between the solution with a very small time step and the analytical solution is regarded as the spatial error, since the error is dominated by the spatial error with a very small time step.

Analysis of the Adaptive Newton Tolerance
This section compares the performance of the proposed Newton tolerance, given by Eq. (37) and = 0.1 , with commonly adopted tolerances given by Eq.  give a reasonable tolerance and still maintain temporal accuracy. In this case, we employ P = 2 and N = 3 with h = 1∕3 and three different time steps: Δt = 0.05, 0.1, and 0.2. Figure 2 shows that large values of the parameters 1 and 2 in the Newton tolerances from Eq. (46) lead to a degradation of the order of temporal accuracy. Small values of these parameters achieve desired temporal order of accuracy but are not necessarily the best choice or the most efficient choice. Furthermore, they need to be calibrated for each  Fig. 2, the adaptive Newton tolerance of Eq. (37) achieves the optimal order of accuracy.
To verify that the Newton tolerance is not unnecessarily small, we numerically seek the maximum Newton tolerance max which maintains a difference between ‖ n+1 t + Δt̄ it ‖ and ‖ n+1 t ‖ is within 1% . Compared with max , the designed adaptive Newton tolerance ‖ t ‖ should be no larger to maintain temporal accuracy and should not be too small to maintain efficiency. Table 1 shows that the ratios between the adaptive Newton tolerance and the max is of the same order of 0.1 for different time steps. This indicates that the iterative error scales simultaneously with temporal error and = 0.1 is close to the maximum choice that can maintain temporal accuracy in this test.

Analysis of Errors in the Adaptive Time Stepping
We described in Sect. 3.1.1 how to chose the tolerances in Newton iteration to ensure that the iterative error is smaller than the temporal error. This section aims at studying spatial, temporal, and total discretization errors characterized in terms of parameters such as the CFL number and the mesh size, h. The adaptive time-stepping strategy is evaluated at a polynomial order P = 6 . The total discretization errors and temporal errors are calculated with respect to the analytical solution and a reference solution computed at a CFL number of 0.01.
The total errors of the third-order ESDIRK (ESDIRK3) with the adaptive time-stepping, CFL = 10.0, CFL = 1.0, and CFL = 0.2 using different mesh sizes, h = 0.2 , h = 0.1 , h = 0.067 and h = 0.05 , are presented in Fig. 3. The corresponding temporal errors are shown in Fig. 4. The spatial errors are also presented for comparison. The temporal errors with the adaptive time steps are always smaller than the spatial errors, as per design, and the total error converges at the same rate as the spatial error. However, for simulations with CFL = 10.0, the total errors are dominated by the temporal errors and converge at a lower order of 3 instead. For simulations with CFL = 0.2, the total errors are dominated by spatial errors. However, the temporal errors are much smaller than the spatial errors, which may lead to extra costs without obviously improving the results. For simulations with CFL = 1.0, temporal error is smaller than the spatial error for h = 0.2 but larger than that for h = 0.05 , which emphasizes the need for calibration for each time step size when using a time step based on a CFL number. Referring to the adaptive time-stepping methods in references [5,22], the temporal error curve should be horizontal lines because the temporal error tolerance in these methods is user defined and independent of h.
The CPU time of the above simulations is compared in Fig. 5. The simulation with the adaptive strategy is close to the most efficient ones at all the error levels. It is slightly more expensive than the most efficient one at some error levels partly due to the extra costs in estimating the errors. Figure 5 also shows the corresponding curves of GMRES iteration number and residual evaluation number. These numbers are the summation of all the time steps during the time interval. The GMRES iteration number and the residual evaluation number of the adaptive method are also close to the smallest among different time-stepping methods. For the parameters adopted in this test, the efficiency of the adaptive stepping can be up to ten times more efficient at the same error level than the simulations with a naive choice of CFL, which highlights the importance of balancing the spatial and temporal errors from the point view of efficiency.

(c) GMRES iteration numbers
The errors of ESDIRK schemes with N = 2, 3, 4 are compared in Fig. 6. All temporal errors are relatively smaller than the spatial errors and the total errors almost coincide with each other for different values of N. Therefore, the simulation results are almost independent of the ESDIRK scheme adopted in the adaptive time-stepping. This property should also be achieved by the adaptive methods proposed in [5,22], but it requires an additional calibrating process which is problem dependent.

Steady-State Flat Plate Boundary-Layer Flow
The objective of this example is to illustrate how the adaptive time-stepping strategy could be used to accelerate convergence towards steady-state in a practical case. The case considered here is the boundary-layer flow past a flat plate at a Mach number Ma = 0.1 and a Reynolds number Re = 1.6 × 10 6 . The geometry and boundary conditions are depicted in Fig. 7. The details of boundary conditions can be found in references [20,29]. The simulations are performed using P = 2 and N = 3 with the highly stretched mesh shown in Fig. 8a. The computed profile of the horizontal velocity is shown and compared with the analytical Blasius profile in Fig. 8b. In all the simulations below, the time step is set to be no larger than that corresponding to a CFL number of 5 × 10 4 to maintain the stability of the simulations.  Eq. (46) is adopted since the default choice of 2 = 10 −3 easily leads to blowing up at such a large CFL number. The growing CFL number method is a commonly adopted strategy for steady-state simulations because a small time step is needed to maintain stability in the initial transient phase, while a large time step is needed to accelerate the convergence when the flow field is close to the steady-state solution [6,28]. The convergence history of the normalized residual norm ‖L ( )‖ is presented in Fig. 9. The simulation with the adaptive time-stepping proposed in Sect. 2.2 converges almost as fast as the simulations with a CFL growth rate of 2.0. The residual, Δt and CPU time at different steps are illustrated in Fig. 10. It shows that the Δt of the adaptive method grows much slower at the first 100 time steps. However, the extra CPU cost because of the extra steps is low since it is cheap to converge at a small Δt . Figure 10 also demonstrates that the adaptive Δt is small at the beginning of the simulation. This is because the transient flow is highly unsteady leading to a large temporal derivative term D t ( ) in Eq. (15) and consequently to a smaller Δt . As the flow field becomes closer to the steady state, the Δt increases because the magnitude of the unsteady waves in the flow field becomes smaller, which leads to a smaller term D t ( ) and larger Δt based on the relation of Eq. (29). The adaptive time-stepping adjusts the Δt based on unsteady properties of the flow field and achieves fast convergence while maintaining stability.
Although the adaptive time-stepping may not be as efficient as methods specially designed for steady-state simulations, such as local time-stepping and multi-grid, it is general and, with a fixed set of parameters, is capable to efficiently simulate a wide range of steady and unsteady flows.

Taylor-Green Vortex
Here, we apply the adaptive time-stepping strategy to the simulation of the Taylor-Green vortex, which represents an unsteady flow of decaying vortices. The initial large vortices break down into smaller scale vortices and eventually into turbulent flow. The simulations are run with Ma = 0.1 , Re = 1 600 and periodic boundary conditions on a 32 3 uniformly distributed mesh with P = 4 and N = 4 . Note that, to simulate nearly incompressible conditions, the flow is taken to be compressible but at a low Mach number of Ma = 0.1. Figure 11 presents the evolution of 2 ∕ , where denotes the enstrophy. According to reference [7], this variable is a good estimate of the kinetic energy dissipation rate in the incompressible limit. Figure 12 presents the errors of the dissipation rate as defined by Eq. (18) of reference [7]. The simulations are run with the set of values = 0.1 and = 0.01 for the adaptive time-stepping strategy and they are also run with Δt = 0.1 and Δt = 0.2 for comparison. The decay curves of 2 ∕ are in relatively good agreement with each other for all the simulations. Their differences can be analyzed from the distributions of the errors in Fig. 12. The error with adaptive time-stepping and = 0.1 is slightly larger than other errors. The error with = 0.01 almost coincides with the error of Δt = 0.1 . Both simulations with = 0.1 and = 0.01 are accurate in temporal accuracy especially compared with the spatial error, which can be roughly estimated by the deviation from the DNS result [7].
We enforce the local temporal error to be smaller than the spatial error by design, but there is no reason to introduce a bias towards the spatial error if the main consideration    is simply to minimize the error magnitude. This consideration may be of relevance for under-resolved turbulent simulations where the spatial truncation error term sometimes serve as an implicit modeling of the under-resolved scales, which is referred to as the implicit large eddy simulation model [9].
To avoid the temporal error from influencing the modeling of the under-resolved scales, a smaller temporal error compared with the spatial error, = 0.01 , is adopted in implicit LES simulations of turbulent flows in our study. This approach gives reasonable results as shown in Fig. 12 and in Sect. 3.4. The computational CPU time of the simulations with = 0.01 is slightly larger than that of simulation with Δt = 0.1 (1.27 times) partly because of the additional costs required for the error estimation.

Turbulent Flow Over a Circular Cylinder at Re = 3 900
The widely studied test case of the turbulent flow over a circular cylinder at Re = 3 900 is adopted here to further illustrate the performance of the adaptive strategy in wall-bounded turbulence simulations. Details of the settings of the test case can be found in reference [29], which employ 123 360 curved unstructured hexahedral elements and a polynomial order of P = 2 . The unstructured mesh distributions and the vortex structures behind the  Fig. 13. Reference [29] adopted a second-order DIRK scheme (DIRK2) with a time step of 0.01, which corresponds to a CFL number of about 40. The time step was determined based on some numerical tests on its efficiency and accuracy in reference [29]. Here, we simulate the same case using the ESDIRK2, ESDIRK3 and ESDIRK4, together with the adaptive time-stepping strategy, to demonstrate their performance.
The adaptive time steps for a duration corresponding to approximately 19 vortex shedding periods are shown in Fig. 14. Table 2 presents the corresponding averages, Δt . The averaged time step of ESDIRK2 is approximately 0.027, which is 2.7 times larger than the Δt = 0.01 of the DIRK2. The adapted time steps are larger for the ESDIRK3 and ESDIRK4 with averaged values of 0.043 and 0.127, respectively, because of the improved resolution of the high-order ESDIRK schemes. Table 2 summarizes the averaged time steps and CPU time of the simulations. The results show that the simulations with adaptive time-stepping strategy achieve good efficiency with only a small overhead in CPU time partly due to the extra costs of error estimations. The time-averaged velocity distributions are presented in Fig. 15. The result of ESDIRK2 with Δt = 0.027 almost coincides with the result of DIRK2 with Δt = 0.01 , both of which are in a good agreement with experimental results [23]. The simulations with adaptive time step give accurate predictions of the time-averaged quantities. Figure 16 depicts the distribution of elemental Δt for the ESDIRK2 scheme corresponding to the variable u on a slice of the flow field. To highlight the elements with a smaller elemental time step than the global time step, the contour coloring is cut below the global time step and these elements are presented in white. One important observation is that the relation of Eq. (29) is maintained in the majority of elements because the elemental time step is smaller than the global time step in fewer than 1% of the elements. The influence of m in Eq. (33) is negligible in the global time step in this test case, which is only effective in elements with very small error magnitudes.

Summary of Parameters in the Simulations
The tunable parameters in the adaptive strategy are listed in Table 3. The parameter controls the relation between the temporal error and spatial error, which is set to 0.1 for general simulations, but to 0.01 for implicit large-eddy simulations as discussed in Sect. 3.3. Similarly, controls the relation between the temporal error and iterative error, which is set to 0.1. The value r in Eq. (35) determines the contribution of elemental time steps to the global time step. For r = 0 , a simple algebraic average is used. A larger r leads to a larger contribution from elements with larger temporal errors. The term m in Eq. (33) avoids division by zero and also serves as a threshold values for the error level of interest. The current choice of Eq. (34) is motivated by the truncation error of the time integration process. In determining and , we assume the spatial error is dominating. However, no optimization on the values of and is performed, since the optimal values are usually problem dependent. On the contrary, in the simulations without the adaptive strategy, the choices of Δt and Newton tolerance have to be changed case by case to achieve a good performance as summarizes in Table 4. These case by case choices are obtained at costs of repeated numerical experiments, which could be expensive especially when starting the simulations of a new test case. The comparison of computational efficiency in Table 4 shows that the adaptive strategy is close to the most efficient choices in all the cases without case by case optimized parameters, which demonstrates the remarkable generality of the adaptive strategy.

Discussion and Conclusions
We have proposed a balanced adaptive time-stepping strategy based on the idea of balancing the errors generated within one time step which has been implemented in an implicit DG high-order spectral/hp element solver.
This adaptive time-stepping strategy maintains temporal accuracy in the sense that the total error is dominated by the spatial error and further decreasing the temporal error will not obviously improve the results of the discrete PDE system, which is verified by numerical experiments in the isentropic vortex, Taylor-Green vortex and turbulent flow over a circular cylinder problem. Moreover, this adaptive time step is a relatively efficient choice because it is around the maximum value with temporal accuracy, the efficiency of which is tested in a variety of steady-state and unsteady problems, including the isentropic vortex, Taylor-Green vortex, flat plate boundary layer and turbulent flow over a circular cylinder. In all these tests, the adaptive methods are close to the most efficient one possibly with small overhead partly due to extra costs in error estimating.
This adaptive time-stepping strategy performs well in a variety of problems without the need of tuning the parameters or requiring a priori knowledge of the flow properties, thus reducing the necessity for user intervention. This feature will facilitate the development of automatic CFD simulation pipelines in a variety of application areas such as, for instance, shape optimization.
The main idea of the adaptive time stepping is to construct a proper relation between different errors, which is not limited to the spatial discretization methods, temporal discretization methods and error estimators adopted in this paper. Provided properly defined spatial and temporal error estimators, the idea should also be applicable to other unsteady PDE solvers. Currently, the adaptive method is based on the local errors generated within one time step, which in theory can not guarantee the global error is still spatial error dominated. This could be improved by adopting some global error estimators.
The adaptive strategy provides an upper bound of the time step based on the requirement of temporal accuracy. However, this upper bound does not necessarily maintain simulation stability in challenging problems. Methods for improving robustness could be considered in the future.

Appendix A Butcher Tableau of ESDIRK Schemes
The coefficients of ESDIRK schemes adopted in this paper can be expressed in the form of the general Butcher tableau shown in Table A1 [17], where S and R = S + 1 are the numbers of stages for the main scheme and the embedded scheme, respectively. The coefficients of the matrix A are sufficient to define the whole Butcher tableau since c i = ∑ S+1 j=1 a ij , b i = a Si , and b i = a Ri . The last two relations indicate that the main and    Tables A2, A3 and A4, respectively.