Spatial coupling of an explicit temporal adaptive integration scheme with an implicit time integration scheme

TL;DR

This paper introduces a spatially coupled hybrid time integration scheme combining explicit and implicit methods, extended to adaptive time stepping on unstructured grids, ensuring stability, accuracy, and conservation.

Contribution

The paper develops an extended hybrid scheme that adaptively couples explicit and implicit time integration methods on unstructured grids with local mesh refinement.

Findings

Linearly stable and second-order accurate scheme.

Validated with wave packet, Sod's tube, and vortex transport tests.

Maintains conservation across interfaces of different cell ranks.

Abstract

The Reynolds-Averaged Navier-Stokes equations and the Large-Eddy Simulation equations can be coupled using a transition function to switch from a set of equations applied in some areas of a domain to the other set in the other part of the domain. Following this idea, different time integration schemes can be coupled. In this context, we developed a hybrid time integration scheme that spatially couples the explicit scheme of Heun and the implicit scheme of Crank and Nicolson using a dedicated transition function. This scheme is linearly stable and second-order accurate. In this paper, an extension of this hybrid scheme is introduced to deal with a temporal adaptive procedure. The idea is to treat the time integration procedure with unstructured grids as it is performed with Cartesian grids with local mesh refinement. Depending on its characteristic size, each mesh cell is assigned a rank. And for two cells from two consecutive ranks, the ratio of the associated time steps for time marching the solutions is $2$. As a consequence, the cells with the lowest rank iterate more than the other ones to reach the same physical time. In a finite-volume context, a key ingredient is to keep the conservation property for the interfaces that separate two cells of different ranks. After introducing the different schemes, the paper recalls briefly the coupling procedure, and details the extension to the temporal adaptive procedure. The new time integration scheme is validated with the propagation of 1D wave packet, the Sod's tube, and the transport of a bi-dimensional vortex in an uniform flow.

Tables

Table 1. Table 1 : Flux formula depending on neighbour cells. The table is ”symmetrical”: the flux between two cells is independent of the direction of information propagation.

	Heun	Hybrid	IRK2
Heun	$F^{Heun}$	$F^{Heun}$	$\times$
Hybrid	$F^{Heun}$	$F^{Hybrid}$	$F^{IRK2}$
IRK2	$\times$	$F^{IRK2}$	$F^{IRK2}$

Table 2. Table 2 : Total error slopes

Time integrator	$E_{tot}$ slope
Heun	1.97
Heun + TA	1.93

Table 3. Table 3 : Total error slopes

Time integrator	$E_{tot}$ slope
Heun	1.97
Heun + TA	1.93
AION	2.08
AION + TA	2.04

Equations

$$\frac{\partial W}{\partial t}+\nabla \cdot F(W)+\nabla \cdot G(W,\nabla W)=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\iiint }_{{\Omega }_j}W\mathrm{d}\Omega =-{\iint }_{A_j}F(W)\cdot \vec{n}\mathrm{d}S-{\iint }_{A_j}G(W,\nabla W)\cdot \vec{n}\mathrm{d}S,$$

$${\overline{W}}_j=\frac{1}{|{\Omega }_j|}{\iiint }_{{\Omega }_j}W\mathrm{d}\Omega .$$

$$\frac{\mathrm{d}{\overline{W}}_j}{\mathrm{d}t}=R({\overline{W}}_j),$$

$$\left\{\begin{aligned} &\frac{\mathop{}\!\mathrm{d}\overline{W}_{j}(t)}{\mathop{}\!\mathrm{d}t}=R\big{(}{\overline{W}}_{j}(t)\big{)},\quad\forall t\in R^{+},\\ &\overline{W}_{j}(0)=W_{j}^{0},\quad\forall j\in\{1,...,N\},\end{aligned}\right.$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor {stage}:}&\displaystyle\widehat{W}=W^{n}+\Delta tR(W^{n}),\\ \vspace*{-3mm}&\\ \hbox{Corrector {stage}:}&\displaystyle W^{n+1}=W^{n}+\frac{\Delta t}{2}\big{(}R(W^{n})+R(\widehat{W})\big{)}.\\ \end{array}\right.$$

$$\displaystyle W^{n+1}=W^{n}+\frac{\Delta t}{2}\big{(}R(W^{n})+R(W^{n+1})\big{)}.$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor stage:}&\displaystyle\widehat{W_{j}}=W_{j}^{n}+\Delta tR(W_{j}^{n})\\ \vspace*{-3mm}&\\ \hbox{Corrector stage:}&\displaystyle\left\{\begin{array}[]{ll}\displaystyle W^{n+1}_{j-1}=W^{n}_{j-1}+\frac{\Delta t}{2}\big{(}F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}}-F_{j-\frac{3}{2}}^{n}-\widehat{F}_{j-\frac{3}{2}}\big{)}\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j}=W^{n}_{j}+\Delta t\big{(}F_{j+\frac{1}{2}}^{Hybrid}-\frac{1}{2}(F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}})\big{)}\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+1}=W^{n}_{j+1}+\Delta t\big{(}\frac{1}{2}(F_{j+\frac{3}{2}}^{n}+F_{j+\frac{3}{2}}^{n+1})-F_{j+\frac{1}{2}}^{Hybrid}\big{)}\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+2}=W^{n}_{j+2}+\frac{\Delta t}{2}(F_{j+\frac{5}{2}}^{n}+F_{j+\frac{5}{2}}^{n+1}-F_{j+\frac{3}{2}}^{n}-F_{j+\frac{3}{2}}^{n+1}).\end{array}\right.\end{array}\right.$$

$$\begin{array}[]{ll}V^{L}=&\omega_{j}\displaystyle\bigg{[}\frac{V^{n}_{j}+\widehat{V}_{j}}{2}\bigg{]}+\frac{1}{2}(\widetilde{\nabla V})_{j}^{n}\cdot\overrightarrow{C_{j}C_{f}}+\bigg{(}\omega_{j}-\frac{1}{2}\bigg{)}(\widetilde{\widehat{\nabla V}})_{j}\cdot\overrightarrow{C_{j}C_{f}}+\\ &(1-\omega_{j})\displaystyle\bigg{[}V^{n+1}_{j}+(\widetilde{\nabla V})_{j}^{n+1}\cdot\overrightarrow{C_{j}C_{f}}-\frac{1-\omega_{j}}{2}(\widetilde{\Delta_{t}V})^{n+1}_{j}\bigg{]},\\ V^{R}=&\omega_{j+1}\displaystyle\bigg{[}\frac{V^{n}_{j+1}+\widehat{V}_{j+1}}{2}\bigg{]}+\frac{1}{2}(\widetilde{\nabla V})_{j+1}^{n}\cdot\overrightarrow{C_{j+1}C_{f}}+\bigg{(}\omega_{j+1}-\frac{1}{2}\bigg{)}(\widetilde{\widehat{\nabla V}})_{j+1}\cdot\overrightarrow{C_{j+1}C_{f}}+\\ &(1-\omega_{j+1})\displaystyle\bigg{[}V^{n+1}_{j+1}+(\widetilde{\nabla V})_{j+1}^{n+1}\cdot\overrightarrow{C_{j+1}C_{f}}-\frac{1-\omega_{i}}{2}(\widetilde{\Delta_{t}V})^{n+1}_{j+1}\bigg{]},\\ \end{array}$$

$$\Delta t=\underset{j\in \mathcal{D}}{\min }\hat{\nu }\frac{h_j}{\Vert \overrightarrow{v_j}\Vert +c_j},$$

$$\displaystyle\omega_{j}=\min\big{(}1,\frac{1}{\nu_{j}}\big{)},\hbox{ with }\nu_{j}=\frac{(\|\vec{v_{j}}\|+c_{j})\,\Delta t}{h_{j}}.$$

$$(\widetilde{\Delta_{t}V})^{n+1,s}_{j}=\left\{\begin{array}[]{ll}\displaystyle\max\left[0,\min\left(\min_{j}[\beta(V_{j}^{n+1,s}-V_{j}^{n})+(\widetilde{\Delta_{t}V})^{n+1,s-1}_{j}],V_{j}^{n+1,s}-V_{j}^{n}\right)\right]&\\ &\hbox{ if }V_{j}^{n+1,s}-V_{j}^{n}\geq 0,\\ &\vspace*{-0mm}\\ \displaystyle\min\left[0,\max\left(\max_{j}[\beta(V_{j}^{n+1,s}-V_{j}^{n})+(\widetilde{\Delta_{t}V})^{n+1,s-1}_{j}],V_{j}^{n+1,s}-V_{j}^{n}\right)\right]&\\ &\hbox{ if }V_{j}^{n+1,s}-V_{j}^{n}<0.\\ \end{array}\right.$$

$$\beta=\frac{2\big{(}1-\omega_{j}^{s}\nu_{j}^{s}(2-\omega_{j}^{s}\nu_{j}^{s})\big{)}}{\nu_{j}^{2}(1-\omega_{i}^{s})^{2}}$$

$$K=\left\lfloor\frac{\ln\Big{(}\frac{\Delta\tau_{j}}{\Delta t_{\min}}\Big{)}}{\ln(2)}\right\rfloor,$$

$$\Delta t_K=2^K\Delta t_{\min }.$$

$${\hat{W}}_{\beta }^{2\Delta t}$$

$${\hat{W}}_{\beta -1}^{2\Delta t}$$

$$F_{\beta+1/2}(\widehat{W}^{\Delta t})=\frac{1}{2}\big{(}F_{\beta+1/2}(W^{0})+F_{\beta+1/2}(\widehat{W}^{2\Delta t})\big{)}.$$

$$R_{\beta +1}^{\Delta t}=\frac{1}{4}{R}_{\beta +1}({\hat{W}}^{2\Delta t})+\frac{3}{4}{R}_{\beta +1}(W^0),$$

$${\hat{W}}_{\beta +1}^{\Delta t}=W_{\beta +1}^0+\Delta t{R}_{\beta +1}^{\Delta t}.$$

$${\hat{W}}_{\beta +1}^{2\Delta t}$$

$$W_{\beta +1}^{2\Delta t}$$

$${\hat{W}}_{\beta }^{\Delta t}$$

$$W_{\beta }^{\Delta t}$$

$${\hat{W}}_{\beta }^{2\Delta t}$$

$$W_{\beta }^{2\Delta t}$$

$${W}^{2\Delta t}_{\beta}={W}^{0}_{\beta}+\frac{\Delta t}{2\Delta x}\bigg{(}2\big{(}F_{\beta+\frac{1}{2}}(\widehat{W}^{2\Delta t})+F_{\beta+\frac{1}{2}}(W^{0})\big{)}-F_{\beta-\frac{1}{2}}(W^{0})-F_{\beta-\frac{1}{2}}(\widehat{W}^{\Delta t})-F_{\beta-\frac{1}{2}}(W^{\Delta t})-F_{\beta-\frac{1}{2}}(\widehat{W}^{2\Delta t})\bigg{)}.\\$$

$$\Delta x\frac{d{\mathcal{W}}_j}{dt}$$

$$\Delta x\frac{d{W}_j}{dt}$$

$$e(W^{2\Delta t})={\mathcal{W}}^{2\Delta t}-W^{2\Delta t},$$

$$F_{j+1/2}$$

$$\frac{{\hat{W}}^{\Delta t}-W^0}{\Delta t}$$

Full text

1 Introduction

The Reynolds-Averaged Navier-Stokes (RANS) equations account for the mean effects of the turbulence on the main conservative quantities. Of course, the RANS method is unable to represent unsteady effects of the turbulence on the flow. But this method has a relatively low CPU cost, and it is accurate enough for the computation of boundary layers for example. For this reason of low CPU cost, it is today one of the preferred technique for use in an industrial context. RANS equations are generally time-marched using an implicit time formulation for fast convergence to the steady-state solution. For unsteady RANS equations, the implicit time integration enables large time steps. Implicit time integration requires to solve large linear system of equations, and then tends to be expensive in terms of CPU cost.

Large Eddy Simulation (LES) solves some of the turbulence effects, and is increasingly being considered in an industrial context. The principle of LES is to introduce a model for the smallest turbulent scales that have a universal nature, and to capture the largest turbulence scales. In practice, the separation between the tractable and the modeled spectra is defined by the mesh spacing and the numerical scheme. Then, LES equations are generally integrated using explicit schemes that exhibit good spectral properties of dissipation and dispersion and can attain any order of accuracy. However, the associated time steps are limited by the Courant-Friedrichs-Lewy (CFL) number.

In [16] was introduced the AION time integration scheme that spatially couples the explicit Heun’s scheme and the implicit Crank-Nicolson’ scheme. The principle of the AION scheme is to blend both schemes using a transition function $\omega$. Grid cells can be declared as explicit, implicit, or hybrid depending on their values of $\omega$, and are associated with their own time integration scheme. Explicit and implicit cells are respectively handled by the Heun’s and Crank-Nicolson’ schemes, while a hybrid scheme is in charge of the other grid cells. An important property of the transition function is to allow a quick switch between the explicit and implicit cells, keeping the number of hybrid cells as low as possible. In addition, attention was paid on the stability analysis to avoid wave amplification for any wavenumber [16]. This paper deals with the adaptation of the AION scheme to a temporal adaptive procedure.

The Adaptive Mesh Refinement (AMR) method allows more grid cells to be placed in regions of interest to better capture the flow physics while keeping a small number of cells in zones of low interest. Many efforts were dedicated to the spatial adaptation, especially for Cartesian grids. For unsteady simulations, it can be of great interest to couple the spatial adaptive method with a temporal adaptive approach. Indeed, a standard unsteady explicit computation is constrained by the CFL condition over the whole computational domain, and the maximum stable time step is generally associated with the smallest cells. It follows that the time integration of the largest cells is performed by a fraction of the maximal local allowable time step, and this leads to useless computations. The principle of temporal adaptive schemes is to allow the mesh cells to integrate the solution using their own time step according to the local CFL condition. This type of local time-stepping method is included in FLUSEPA11footnotemark: 1 solver [6]. This solver, developed by ArianeGroup, is used in the certification process of launchers, when complex physical phenomena occur.

Concerning AMR, the main difficulty lies in the treatment of flux at the interface between cells of different grid sizes. Sanders and Osher [17] proposed a one-dimensional local time-stepping algorithm and later on, Dawson [9] extended the procedure to multi-dimensional problems. The procedure is built in order to keep a space time conservation: for an interface between large and small cells, the sum of fluxes on small cells is exactly the one on the larger cell. However, the procedure leads to a first-order time scheme. An extension to the second order of accuracy was proposed more recently. Dawson and Kirby [10] extended the procedure by means of Total Variation Diminishing (TVD) Runge-Kutta schemes. Constantinescu and Sandu [7] developed a set of multirate Runge-Kutta integrators called Partitioned Runge-Kutta. They are adapted to automatic mesh refinement, and respect strong mathematical properties (Strong Stability Preserving). Later on, the same authors extended the procedure to the multirate explicit Adams scheme [21].

Berger and Oliger [4] and Berger and Colella [3] proposed a different approach for playing with both space and time refinements. Starting from the initial grid level indexed $l=0$, several grid levels $l=1$,…, $l=l_{max}$ are introduced by local refinement (see Fig. 1). For transferring information from coarse to fine grids, ghost cells are introduced at the level $l$, and the values inside them are interpolated from the coarse grid at level $l-1$. It must be highlighted that these ghost cells can be found on boundary conditions or inside the computational domain, and serve as local boundary conditions for the grid level. The interest of such method lies in its simple implementation.

The procedure consists of time integrating the mesh by means of the time classes, starting from a coarse grid ($l=0$) and finishing by the most refined level ($l=l_{max}$). On any grid level, the following steps are considered:

• •

time integrate cells of level $l$ with local time step ($\Delta t_{l}$). If $l>0$, the boundary information is collected from the next coarser level $l-1$ by interpolation of the flux.

• •

synchronise cells of levels $l$ and $l+1$, and interpolate correction to finest cells (of level $[l+2,...,l_{max}]$).

Hence the adaptive scheme allows the time integrations of the different level of grids by means of their own time steps, and more time integration steps are needed on the refined levels than on the coarser one. The synchronization step must be seen as a correction that makes the procedure consistent. Nevertheless, the procedure must be defined carefully since synchronization and interpolation drive order of accuracy and conservation. For instance, let’s take the example introduced in Fig. 2. The principle of “internal” ghost-cell $C_{j_{\alpha}}$ is to serve as an additional variable. The ghost cell value and the gradient can be interpolated from cell-centered values in the surrounding cells. However, this can lead to a non conservative procedure since the interface state computed from the ghost cells (using state and gradients for instance) and the one from the large cell $C_{j}$ may be different.

Several methods were introduced to limit the impact of loss of accuracy. Bell et al. [2] introduced a correction on the state $W$ thanks to a passive scalar, which allows a fast convergence for quasi-steady state. But they noticed that the conservation of flux is not guaranteed. Bell [1] also applied a correction seen as a kind of ”fixed” Dirichlet boundary condition but involving an inconsistency at interfaces between grids.

The temporal adaptive scheme in FLUSEPA follows some basic principles introduced by Kleb et al. [13]. Their procedure is based on the explicit Euler time integration, and leads to a first-order time-accurate solution. Obtaining a second-order time integration is a prerequisite, and the approach shares some ideas with Krivodonova’s work [14], even if the procedure differs in the way the time integration is performed on classes holding the smallest cells. The first specific point is the use of the second order explicit scheme proposed by Heun. The second point of importance is the extension to the adaptive time treatment for Heun’s scheme. The last key ingredient is the coupling of the adaptive time integration with the hybrid AION scheme presented in [16]. The remainder of this paper unfolds as follows. In Sec. 2, the standard form of the Navier-Stokes is recalled and the explicit, implicit and hybrid schemes chosen for the present study are introduced. Our explicit, implicit and hybrid basic schemes (Heun’s, Crank-Nicolson’s and AION schemes) are defined in Sec. 3. In Sec. 4, the temporal adaptive method implemented in FLUSEPA is introduced. It must be underlined that the adaptive version of the Heun’s scheme was briefly described in [5], and the full description of the scheme implemented in FLUSEPA is also one objective of the paper. Mathematical analysis including order of accuracy, stability and spectral behaviour is then presented. In Sec. 6, the procedure is extended to our coupled explicit / implicit time integrator called AION. Before concluding, Sec. 8 is dedicated to the validation of the adaptive AION scheme for 1D and 2D test cases with temporal adaptive approach.

2 Discretization of the Navier-Stokes Equations

The Navier-Stokes system of equations is written in the following compact conservation form,

[TABLE]

with $W$ the vector of conservative variables, its gradient $\nabla W$, the convection flux $F$ and the diffusion flux $G$. The $N$ non-overlapping rigid stationary cells $\Omega_{j}$ map the computational domain $\Omega$, and Eq. (1) is integrated over every mesh cell. Applying the Gauss relation that ties the volume integrals of the divergence terms to the interface fluxes gives the weak form,

[TABLE]

where $\Omega_{j}$ is the $j$-th control volume with its border $A_{j}$, and $\vec{n}$ is the outgoing unit normal vector. The averaged conservative variables are defined as

[TABLE]

This relation (3) allows to rewrite the conservation laws (2) discretized with a finite-volume formulation in the following differential form,

[TABLE]

where $R(\overline{W}_{j})$ is the residual computed using the averaged quantities $\overline{W}_{j}$ in cell $j$. For any cell, the residual is the sum of the flux over the whole boundary of a cell.

In this paper, convection and diffusion fluxes are discretized by means of a $k$-exact formulation coupled with successive corrections [19, 11]. The principle of the formulation is to define a local polynomial approximation of the unknowns. The successive corrections are designed to avoid geometrical reconstructions for parallel computations. Starting from the solution, the Taylor expansion of the unknowns defines the local polynomial reconstruction, and the order of accuracy is defined by the error term in the Taylor approximation. The process for defining the coefficients of the Taylor expansion can be found in [19, 11].

Finally, the finite-volume formulation of Eq. (4) can be expressed as the following Cauchy problem,

[TABLE]

where ${W_{j}^{0}}_{1\leq j\leq N}$ denotes the initial solution in the mesh cells. For the sake of clarity, the averaging symbol will be dropped, and ${W}$ will represent the averaged quantities over the control volumes.

3 Heun, Crank-Nicolson, and AION methods for time integration

Heun’s, Crank-Nicolson’s and AION schemes, which are considered as standard time integration schemes in this study, are presented below.

3.1 Heun’s scheme

Heun’s explicit scheme [12] is a second-order accurate predictor-corrector method. The state $W^{n+1}$ is first predicted with a forward Euler scheme, and then corrected by a standard trapezoidal rule,

[TABLE]

3.2 Crank-Nicolson’s scheme (IRK2)

Crank-Nicolson’s implicit scheme [8] is a second-order accurate method defined as

[TABLE]

Eq. (7) is solved with an iterative Newton method since the residual $R$ is nonlinear in $W$. This scheme belongs to the class of second-order Implicit Runge-Kutta schemes, and in the following, it will also be denoted IRK2 for conciseness.

3.3 The AION scheme

The AION scheme couples Heun’s and Crank-Nicolson’s schemes by means of a hybrid scheme that enables a smooth transition between the latter two schemes. The AION scheme is defined to ensure a unique flux on each interface, leading to a conservative formulation. The full scheme description is available in [16] and the key points are recalled below for the sake of clarity.

The cell status $\omega_{j}$ is used to smoothly switch from a time-explicit ($\omega_{j}=1$) to a time-hybrid ($0.6<\omega_{j}<1$) and then to a time-implicit ($\omega_{j}\leq 0.6$) scheme. The cell status is named Heun for $\omega_{j}=1$, IRK2 for $\omega_{j}\leq 0.6$, and hybrid otherwise.

The flux formula used for an interface separating two cells depends on the status of these cells as shown in Tab. 1. The flux formula type is determined using the cell status provided on Fig. 3.

According to Tab. 1, for the 1D example on Fig. 3, the AION scheme is expressed as:

[TABLE]

The reconstructed primitive states (vector $V$) at an interface between cells with indices $j$ and $j+1$, which are used in the computation of the hybrid flux $F^{Hybrid}$, are defined as

[TABLE]

where $C_{j}$, $C_{j+1}$ and $C_{f}$ represent the centers of cell $j$ and $j+1$ and the interface center, respectively.

Remark: A situation not adressed in Tab. 1 nor in Fig. 3 must be described. For any hybrid cell (according to the value $\omega_{j}$) with hybrid flux for all the faces, the predictor state $\widehat{W}_{j}$ is time-integrated by accounting for $\omega_{j}$: $\widehat{W}_{j}=W_{j}^{n}+\omega_{j}{\Delta t}R(W_{j}^{n})$.

The cell status $\omega_{j}$ is locally adapted to the flow and the stability constraint. First, the global time step is chosen to be the maximum allowable time step on a part $\cal{D}$ of the whole computational domain without source of stiff terms,

[TABLE]

where $\hat{\nu}<1$, $h_{j}$ is a reference length scale, $\vec{v_{j}}$ the velocity vector, and $c_{j}$ the speed of sound in cell $j$. Then, the parameter $\omega_{j}$ is defined as

[TABLE]

In Eq. (9), the term ${(\widetilde{\Delta_{t}V})}^{n+1}_{j}$ comes from a generalization of the TVD property to the coupled space/time reconstruction. All details regarding this specific term and the proof of the TVD property are provided in previous paper  [16]. It can be seen as a time slope limiter. This term is given by a Newton’s algorithm which at any step $s$ reads,

[TABLE]

The parameter $\beta$ is computed from the local CFL number $\nu_{j}$ and $\omega_{j}$ at each step $s$,

[TABLE]

This method is second-order accurate in both space and time in fully explicit, fully implicit and explicit-implicit domains.

The next section is dedicated to the time-adaptive extension of the Heun’ scheme.

4 Temporal Adaptive Heun’s scheme

The time-adaptive extension of Heun’s scheme is fully described below. It starts from the very brief introduction in [5].

4.1 Time class for the mesh cells

The local time step of each cell $j$ is computed, $\Delta\tau_{j}=\mbox{CFL}\frac{h_{j}}{\|\vec{v_{j}}\|+c_{j}}$ where $h_{j}$ is a reference length scale, $\vec{v_{j}}$ the velocity vector and $c_{j}$ the speed of sound in cell $j$. The minimum value of the local time step, $\Delta t_{\min}$, enables to define the time class of rank $K$ with

[TABLE]

from which the time step of the class denoted $K$ is deduced,

[TABLE]

From Eq. (15), it is clear that the class of rank [math] is associated with the time step $\Delta t_{\min}$. If the solution of the partial differential equation is integrated in time in the ascending order of class rank, the associated time steps are $2\Delta t_{\min}$, $4\Delta t_{\min}$, and so on. In order to attain the same physical time, cells of class $K$ must be integrated one more time than cells of class $K+1$. In practice, it is also required that the difference in class rank between two cells sharing the same interface is not greater than one. Therefore, two adjacent cells are integrated with the same time step if they are in the same class, or with different time steps with a ratio of 2. There are other approaches with temporal adaptive methodology, where each cell is time integrated with its own maximal allowable time step [25, 24, 22].

In order to obtain time synchronization, attention is paid on the treatment of the flux between cells that belong to classes with different ranks since it is the key point to keep a conservative treatment [4, 3].

4.2 Update of the Solution

When two cells sharing a common interface belong to two different time classes, the cell with the lower rank is iterated one more time than cells with the higher rank. It is then necessary to elaborate a specific computation of the flux on the interface between these two cells to perform time synchronization with conservation, while maintaining the local scheme accuracy.

The method is illustrated on a one-dimensional example with cells of class $K=0$ and $K=1$ in Fig. 4 with space as the abscissa and time as the ordinate. The dashed lines represent the time at which the solution must be computed starting from the initial solution represented by symbols. The first-order time-accurate solutions using the predictor stage of Heun’s scheme will be represented by $\widehat{\times}$.

Starting from the initial time $t=0$, the final time reached by all cells is equal to $2\Delta t$, which is the time step of the class with rank $K=1$ and twice the time step of the cells with rank $K=0$. In the following, $F_{i}^{k\Delta t}$ will represent the flux at the interface $i$ at the time $k\Delta t$. If $k=0$, the superscript will be simply [math]. In the standard finite volume approximation, the local residual of a cell of index $j$ is simply $R_{j}=F_{j+1/2}-F_{j-1/2}$.

Heun’s scheme is an explicit second-order time integrator. Since it contains a predictor-corrector procedure [12], cells with class rank $1$ will undergo only two stages to reach $t=2\Delta t$:

• a-1.

$\displaystyle\widehat{W}^{2\Delta t}=W^{0}+2\Delta tR(W^{0})$

• b-1.

\displaystyle W^{2\Delta t}=W^{0}+\frac{2\Delta t}{2}\big{(}R(W^{0})+R(\widehat{W}^{2\Delta t})\big{)},

but cells with class rank [math] will undergo four stages:

• a-0.

$\displaystyle\widehat{W}^{\Delta t}=W^{0}+\Delta tR(W^{0})$

• b-0.

\displaystyle W^{\Delta t}=W^{0}+\frac{\Delta t}{2}\big{(}R(W^{0})+R(\widehat{W}^{\Delta t})\big{)}

• c-0.

$\displaystyle\widehat{W}^{2\Delta t}=W^{\Delta t}+\Delta tR(W^{\Delta t})$

• d-0.

\displaystyle W^{2\Delta t}=W^{\Delta t}+\frac{\Delta t}{2}\big{(}R(W^{\Delta t})+R(\widehat{W}^{2\Delta t})\big{)}.

$\bullet$ **Step 1:

**The residual $R(W^{0})$ is computed using the initial solution, and the predicted states are obtained for all cells (stages a-0 and a-1), as shown in Fig 5. At this time, it should be highlighted that the stage b-1 needs the residual using $\widehat{W}^{2\Delta t}$ available only in the cells of class 1 (and not in class 0), and the stage b-0 needs the residual using $\widehat{W}^{\Delta t}$ available in the cells of class 0 (and not in class 1).

$\bullet$ **Step 2:

**The residual $R(\widehat{W}^{2\Delta t})$ needed by the stage b-1 needs the predicted states on cells of rank 0 on the required stencil to compute the interface flux (Fig. 6). For the required cells with class rank 0, the estimated states are

[TABLE]

These states are called extrapolated states since $2\Delta t$ violates the CFL stability condition for these cells ($2\Delta t$ is the time step of cells from class of rank $1$).

For cells of class rank $1$, once the flux $F_{\beta+1/2}(\widehat{W}^{2\Delta t})$ is computed, the final state of stage b-1 can be computed at $t=2\Delta t$ (Fig. 7).

$\bullet$ **Step 3:

**For the cells of class rank [math] sharing a face with a cell of class rank $1$, the update of the solution associated with stage b-0 needs the definition of the flux at the interface $\beta+1/2$,

[TABLE]

This interpolation is represented by mean of arrows in Fig. 8.

$\bullet$ **Step 4:

**The flux is directly given by Step 3 for the interface $\beta+1/2$ but the residual in cells $\beta$, $\beta-1$, etc. may involve states in the cells of class 1. The number of cells to treat in class rank 1 depends only on the stencil of the spatial scheme for the cells in class rank 0 (MUSCL formulation for instance). In order to compute the other corrector flux at $t=\Delta t$ for cells of class rank 0 located in the stencil of the spatial schemes (cell $\beta-1$ for instance), it is necessary to predict the states in some cells of class rank $1$ at $t=\Delta t$. A simple first-order accurate prediction would be $\widehat{W}^{\Delta t}_{\beta+1}=W^{0}_{\beta+1}+\Delta tR(W^{0})$. To keep a second-order time accuracy, the residual at time $\Delta t$ is computed with a parabolic interpolation of the residuals,

[TABLE]

which is used to compute the predictor state at time $\Delta t$ using

[TABLE]

The new estimated solution in the cell $\beta+1$ is represented as a predicted state in Fig. 9.

Using the new available information, the cells of class rank 0 can be updated and the stage b-0 is completed. From now on, the variables $W^{\Delta t}$ are available for any cell of class rank 0 (Fig. 10). The last step is to apply stages c-0 and d-0.

$\bullet$ **Step 5:

**Finally, the cells of class rank [math] have to be integrated in time from $\Delta t$ to $2\Delta t$ (stages c-0 and d-0). There are two points to perform the time integration. First, it is mandatory to keep the interface flux constant for the Heun’s stages: the flux computed using Eq. (17) is imposed to compute the residual in cell $\beta$ for stages c-0 and d-0. Moreover, the high-order spatial interpolation (for instance for interface $\beta-1/2$) may need the fields in the cells of class rank 1 at time $\Delta t$. In that case, the data computed using Eq. (19) is kept constant for stages c-0 and d-0.

To conclude, in one time iteration, this temporal adaptive method allows to integrate small cells in time with the time step of the biggest cells of the domain thanks to a subcycling process. Brenner [5] showed that for a maximal time step $\Delta t_{max}$ equal to $2^{K}\Delta t_{min}$, if most cells are included in the time class of rank $K$, then the computational cost is divided by $2^{K}$ in the most favorable case.

4.3 Conservation Property

As stated in Sec. 1, conservation associated to the flux balance is not always guaranteed, and this section is devoted to the proof of conservation for the proposed time-adaptive Heun’s scheme.

First, according to the one-dimensional configuration presented in Sec. 4.2, the time integration of cell $\beta+1$ of class rank $1$ may be formulated with predictor and corrector stages going from $t=0$ to $t=2\Delta t$ as

[TABLE]

All negative contributions from the interface $\beta+\frac{1}{2}$ must be recovered in the flux balance of the cell $\beta$.

Moreover, focusing on the opposite side of the interface, the time integration from $t=0$ to $t=2\Delta t$ of the cell $\beta$ of class rank [math] may be formulated as

[TABLE]

in a first phase, and

[TABLE]

in a second phase. If ${W}^{\Delta t}_{\beta}$ of Eq. (21) is replaced in the second equation of (22), and reminding that Eq. (17) holds, then

[TABLE]

It appears that the time integration of state $W$ at cells $\beta$ and $\beta+1$ between $t=0$ and $t=2\Delta t$ share exactly the same flux at the interface $\beta+\frac{1}{2}$.

5 Accuracy and spectral analysis of temporal adaptive method

5.1 Time Accuracy

In this section, a local analysis of the accuracy for the temporal adaptive version of Heun’s scheme is performed near an interface concerned by time synchronization, which means that the interface is shared by two cells belonging to to adjacent class ranks. For the sake of clarity, the proof is performed in 1D but the demonstration can be kept in a multi-dimensional framework.

Let’s introduce a 1D mesh with a fixed grid size $\Delta x$ with two temporal classes imposed. A generic partial differential equation is integrated spatially, as for standard finite volume. Indeed, two relations are obtained, the first one being associated to the exact relation and the second one to the approximated discrete version:

[TABLE]

with $\mathcal{W}$ the exact state, $\mathcal{R}$ the exact residual and $\mathcal{F}$ the exact flux. As introduced above, $W$ represents the numerical state, $R$ the numerical residual and $F$ the numerical flux. In the following, the numerical error on the state for cells of rank [math] and $1$ between $t=0$ and $t=2\Delta t$ is defined as:

[TABLE]

This numerical error will be computed at each step of the temporal adaptive method presented in Sec. 4.2. From now on, it is assumed that the computation of the flux is $p^{\mbox{\scriptsize th}}$-order accurate in space and that a first-order finite difference approximates the time derivative, which leads to:

[TABLE]

Using Eq. (24), it is clear that the residual behaves as $\mathcal{O}(\Delta x^{p-1})$ due to the relation between the residual and the flux. According to Eq. (26), Eq. (24) can be written as (omitting subscript $j$)

[TABLE]

Equation (27) is typically what occurs during a predictor stage of the Heun’ scheme.

Let us now study the local error performed at cells of class rank $1$.

5.1.1 Local numerical error of cells of class rank $1$

According to the numerical error performed at the predictor state from Eq.(27), and knowing the orders of accuracy of the numerical scheme, the accuracy of the numerical flux (omitting subscript of interface index) at cells of temporal class rank $1$ is:

[TABLE]

Then according to Eqs. (28) and (24), it is possible to formulate the residual as

[TABLE]

With the assumption that no numerical error occurred at instant $t=0$, i.e. $W^{0}=\mathcal{W}^{0}$ and $\mathcal{R}=\mathcal{R}(W^{0})$, the Taylor expansion of flux (omitting subscript of interface index) around ${W}^{0}$ gives:

[TABLE]

The residual is reformulated thanks to Eqs. (30) and (24) as

[TABLE]

Then the state $W$ at cells of class rank $1$ of the corrector step is given by

[TABLE]

Hence it is possible to formulate the numerical error of $W^{2\Delta t}$ at cells of class rank $1$ (see Eq. (25)) using the Taylor expansion of ${\mathcal{W}}^{2\Delta t}$ and Eqs. (31) and (32) as

[TABLE]

5.1.2 Local numerical error due to time interpolation for cells of class rank $1$

The flux at interface $\beta-\frac{1}{2}$ needs a special treatment due to the use of a MUSCL formulation: a special treatment of cells of class rank 1 is performed according to step 4. Here only cell $\beta+1$ of class rank $1$ is considered, so it is important to evaluate the error generated during the parabolic interpolation. Thanks to Eq. (31), the state $\widehat{W}^{\Delta t}_{\beta+1}$ is formulated as:

[TABLE]

Then, the error $e(\widehat{W}^{\Delta t}_{\beta+1})$ according to parabolic interpolation is:

[TABLE]

According to Eq. (35), the numerical flux at interface $\beta-\frac{1}{2}$ is expressed as

[TABLE]

Now let us consider the case of cell $\beta$ of class rank [math].

5.1.3 Local numerical error of cells $\beta$ of class rank [math]

According to the numerical error obtained in Eq. (27) at the predictor step on $\widehat{W}^{2\Delta t}$ and the Taylor expansion of flux around $W$,

[TABLE]

then the flux $F(\widehat{W}^{\Delta t})_{\beta+\frac{1}{2}}$ is written as

[TABLE]

Thanks to Eqs. (36), (38), and (37), the state $W^{\Delta t}_{\beta}$ is formulated as

[TABLE]

And the error $e(W^{\Delta t}_{\beta})$ is written as

[TABLE]

Considering the local error of $W$ at cell $\beta$ until $t=\Delta t$ in Eq. (40), and the error associated to the stencil of the MUSCL reconstruction technique in Eq. (35), the fluxes $F_{\beta-\frac{1}{2}}(W^{\Delta t})$ and $F_{\beta-\frac{1}{2}}(\widehat{W}^{\Delta t})$ are expressed as

[TABLE]

Thanks to Eqs. (41), (38), (37), and (40), the state $W^{2\Delta t}_{\beta}$ is then formulated as

[TABLE]

Then the local error $e({W}^{2\Delta t}_{\beta})$ is

[TABLE]

5.1.4 Partial conclusion

In the case of hyperbolic equations, the CFL condition shows that the time step $\Delta t$ depends linearly on the grid spacing $\Delta x$. Then the local numerical error performed until $t=2\Delta t$ on cell $\beta$ of class rank [math] is $O(\Delta t^{2},\Delta x^{p})$, and the local numerical error at cells of class rank $1$ is $O(\Delta t^{2},\Delta x^{p})$. It appears that the temporal order of accuracy is kept constant locally.

In the following section, the influence of local numerical errors due to the temporal adaptive method on the global error is investigated.

5.2 Numerical assessment of the theoretical behavior

A one-dimensional numerical analysis of the global order of accuracy is performed by solving the advection equation of a sinus wave inside a periodic domain of length $L=1$ ($x\in[0,1]$). The initial state is

[TABLE]

The grid contains $N$ cells, and then the grid size is constant. Two temporal classes are imposed with $100$ cells in the class of rank [math] in the middle of the domain ($\frac{N}{2}-50\leq j\leq\frac{N}{2}+50$). The computation is performed during $3s$ of physical time with a fixed time step $\Delta t=2.5\,10^{-5}s$, with a $1$-exact spatial scheme. The accuracy of the temporal adaptive method with Heun’s scheme (defined as “Heun+TA”) is compared with the standard Heun’s integrator. As expected, the computation of the total errors $E_{tot}$ reveals that the temporal adaptive method conserves the second order space-time accuracy of Heun’s scheme, according to Tab. 2.

5.3 Space-Time von Neumann Analysis

The one-dimensional linear advection equation is considered again in this section:

[TABLE]

with a harmonic initial condition and with periodic boundary conditions,

[TABLE]

with $i^{2}=-1$ and $c>0$. $k$ is the wavenumber of the initial solution. The problem admits an exact solution, and the comparison between the exact theoretical solution and the numerical approximation gives information on both dissipation and dispersion of the fully discrete scheme. The fully discrete relation obtained from Eq. (45) using a second-order finite volume scheme can be written in the following form:

[TABLE]

where $y_{j}^{n}$ represents the harmonic (averaged) solution at discrete time $n$ and in the cell $j$. The complex coefficient $G_{j}$ represents the transfer function between two consecutive time solutions. The coefficient $G_{j}$ can be expressed as $G_{j}=|G_{j}|\,\exp(i\arg(G_{j}))$. In the following, $\mu_{j}=|G_{j}|$ is the dissipation coefficient, and $\phi_{j}=arg(G_{j})/\hbox{CFL}$ the dispersion coefficient. Of course, the definition of $G_{j}$ in the cell $j$ depends on the CFL value, on the class rank and on the spatial accuracy used for the high order definition of the flux (linear extrapolation for second-order scheme). For the sake of clarity, the dependency on all the variables will be omitted.

As presented in the previous section, the temporal adaptive method involves sub-cycling of time integrations. Thus the space-time spectral analysis needs to be performed on several time steps. In the following, it is performed on two time steps for two classes of cells. With $c>0$, two configurations of time synchronisation between classes will be studied as ”step DOWN” and ”step UP” (Fig. 12). For this theoretical analysis the cell size will be fixed to $\Delta x$ in the entire numerical domain. In case of irregular mesh a space-time spectral analysis may be non-trivial [26]. In addition, a numerical analysis of wave packet propagation will be presented in the test case section 8 in order to study the effect of an irregular grid for our time integration setting.

Remark: For a time step $\Delta t$ corresponding to $\mbox{CFL}=0.5$ in cells of class rank [math], the time step $2\Delta t$ corresponds to $\mbox{CFL}=1$ in cells of class rank $1$ since $\Delta x$ is constant in the whole domain. Then the time step $\Delta t$ must be chosen so that $\mbox{CFL}<0.5$ in order to remain below the stability condition in class rank $1$.

In the following spectral analysis, the domain configuration shown in Fig. 12 is considered with $N=300$ cells.

5.3.1 Transfer function for Heun’s scheme, far from class transition

The Heun time integration at cell $\beta-10$ (far from class transition) from $t=0$ to $t=2\Delta t$ is performed such that

[TABLE]

and then the following finite volume formulation is obtained

[TABLE]

5.3.2 Analysis for the cell $\beta$ near step UP

The finite volume formulation of the state $y^{2\Delta t}_{\beta}$ at the cell $\beta$ near the step UP is different:

[TABLE]

Considering the MUSCL reconstruction for a simple advection from left to right in one dimension, the linear extrapolation of the unknown used for the flux computation reads:

[TABLE]

Then, by substitution in Eq. (50),

[TABLE]

Here the terms $\widehat{y}_{\beta}^{2\Delta t}$ and $(\widehat{\nabla y})^{2\Delta t}_{\beta}$, may be source of instability because of the violation of CFL condition for these cells.

5.3.3 Analysis for the cell $\gamma$ near step DOWN

In the configuration of the step DOWN, the time integrations of the state $y$ at the cell $\gamma$ are formulated as:

[TABLE]

5.3.4 Conclusion on the analysis

To obtain the amplification factor $G$ between $y^{2\Delta t}$ and $y^{0}$, all the terms in Eqs. (49), (52), and (53) need to be expressed as a transfer function from $y^{0}$:

[TABLE]

For the step UP configuration, Figs. 13 and 14 show that amplification is never created. For the step DOWN configuration, it appears that the terms $\widehat{y}_{\gamma-1}^{2\Delta t}$ and $(\widehat{\nabla y})^{2\Delta t}_{\gamma-1}$ can create amplification as shown in Figs. 13 and 14.

At CFL$=0.6$, Fig. 15 shows that both step UP and DOWN configurations present strong amplification due to the fixed grid size $\Delta x$ (cf. Remark of Sec. 5.3). It is important to remember that amplification in a certain range of wave numbers does not necessarily imply instability. Nevertheless, this type of configuration is not allowed to appear in a large part of the computational domain. In particular, if several classes are located closely, a wave can be amplified at any step DOWN: CFL below 0.6 is for sure a good solution. Concerning dispersion behaviour, Figs. 13 and  15 reveal that the time synchronisation does no affect the numerical speed of the advected waves for low wave numbers: the dispersion of Heun’s scheme is recovered. However, taking again Heun’s scheme as a reference, Fig. 14 shows that dispersion is changed for wavenumbers above $\frac{\pi}{2\Delta x}$.

5.4 Analysis of $q-$waves

For an one-dimensional advection equation at constant speed $c$, any information is transported at the velocity $c$ and this result is at the core of the theory for (linear) hyperbolic equations. Numerically, numerical parasite waves propagating in the opposite direction of $c$ can be encountered and they are referred as $q-$waves. In addition, a wave propagating in the same direction as $c$ is called $p-$wave. The existence of these $q-$waves was first highlighted by Vichnevetsky [27], and then by many other authors, such as Poinsot and Veynante for combustion [18] or Trefethen [15]. The existence of $q-$waves is not related to the spectral analysis of a numerical scheme. Existence of $q-$waves is linked with numerical group velocity according to Sengupta et al. [23]. In this context, new quantities related to Eq. (45) are introduced according to the definition proposed in [23]. The analysis of $q$-waves for the temporal adaptive method is performed with the time step $2\Delta t$ according to the configuration presented in Sec. 5.3. The numerical phase speed is defined as

[TABLE]

from which is deduced the numerical group velocity,

[TABLE]

For an advection equation with a positive advection velocity $c>0$ and information going downstream, $q-$waves propagate upstream with a negative group velocity $V^{gN}$ [23]. The negative group velocity corresponds to a necessary condition to make $q-$waves appear and propagate, but it is not a sufficient condition. Moreover, the effect of $q-$waves also depends on the amplification factor of the space-time discretization. Indeed with an excessive dissipation or with filtering, the $q$-waves can be damped efficiently. In addition, initial condition, grid resolution and multi-dimensional case have effects on $q$-waves [23]. Here, the analysis is restricted to the one-dimensional case.

Both configurations (steps UP and DOWN) of temporal adaptation coupled with Heun’s time integration are taken into account, as in Sec. 5.3. The numerical phase speed $c^{N}$ and the numerical group velocity $V^{gN}$ are computed from $G$ (Eq. (47)) for cells $\beta$ and $\gamma$. The analysis of $q-$waves is performed for CFL $\in[0,1]$. In Figs. 16, dashed curves correspond to $|V^{gN}|/|c|>>1$, which represents discontinuity of dispersion $\phi$, and grey zones correspond to negative group velocity, $V^{gN}<0$.

Remark: The CFL domain $[0,1]$ is larger than the stable CFL domain, and amplification always occurs for CFL$>0.6$. However, the analysis can still be performed for any value of the CFL number.

Steps UP and DOWN lead to equivalent areas of negative-$V^{gN}$ group velocity (Fig. 16) and the shape of iso contours are also almost the same for CFL $\in[0,0.8]$. In the step UP case (cell $\beta$), another negative-$V^{gN}$ group velocity zone is present for CFL$>0.9$.

The last question concerns the capability of the temporal adaptive Heun scheme to damp $q-$waves. To answer it, it is mandatory to couple the observations on negative group velocity with the dissipation property of the time integration. As mentioned in the remark of Sec.5.3, the numerical scheme is unstable for CFL$>0.5$ and the negative-$V^{gN}$ zone found for CFL$>0.5$ (Fig. 17) matches with the dissipation factor $\mu_{\beta}>3$, leading to amplification. Nevertheless, the zones of negative-$V^{gN}$ for CFL$<0.5$ corresponds to $\mu<0.5$ and the wave is partially dissipated in one iteration.

5.5 Partial conclusion

In this section, the numerical analysis of the proposed time-adaptive version of Heun’s scheme was analysed, focusing on dissipation, dispersion, accuracy and the possibility to encounter $q-$waves.

The next section is dedicated to the very same analysis for our AION scheme that blends the scheme of Crank-Nicolson with the second-order implicit Runge-Kutta integrator.

6 Temporal adaptive method with AION scheme

The standard AION time integration coupled with a 1-exact spatial scheme was shown to be second-order accurate both in space and time [16]. The AION time integrator is built using three underlying schemes applied either in the explicit, implicit or hybrid regions of the mesh. The time step is de facto limited by the CFL condition of the smallest explicit cells. For explicit cells, this time integrator was also designed to match with the predictor-corrector formulation of the standard Heun’s scheme. Indeed, it was highlighted in Sec. 4 that the predictor-corrector formulation is useful for time synchronisation between temporal classes. Hence, the AION time integration may be coupled with the standard temporal adaptive approach. In the latter case, the implicit and hybrid cells may belong to the temporal class of highest rank with the time step limited by the CFL condition of the largest explicit cells. In some cases, it is possible to obtain configurations for which hybrid cells belonging to the highest rank temporal class share an interface with explicit cells of lowest rank temporal classes (even though the approach forbids neighbour cells with a difference in temporal class rank greater than $1$). Thus it is necessary to extend the time adaptive approach to the hybrid scheme of AION time integrator.

In this section, the goal is to propose a way to couple the AION scheme with time adaptation. Time adaptation should not be considered for the treatment of implicit cells since they are protected by the design of the transition parameter. However, it will be mandatory to allow the time synchronisation of temporal classes in hybrid and/or explicit cells. For $\omega=1$, the AION scheme reverts to Heun’s scheme and the very same time synchronisation can be chosen. The situation differs when synchronization occurs in the hybrid part of AION scheme, where only one residual is computed (Eq. 8). In the following, the scheme for which time adaptation and synchronization occur in the hybrid part of the AION scheme will be denoted AION+TA, while Heun+TA will indicate a time adaptation procedure with synchronization in the explicit cells, as presented in Sec. 4.

6.1 Update of the Solution for hybrid time integration

Let us consider again the one-dimensional example introduced in Sec. 4.2. As before, cells with class rank $1$ will undergo only two stages to reach $t=2\Delta t$ thanks to the hybrid time integration:

• a-1.

$\displaystyle\widehat{W}^{2\omega\Delta t}=W^{0}+2\omega\Delta tR(W^{0})$

• b-1.

$\displaystyle W^{2\Delta t}=W^{0}+{2\Delta t}R^{Hybrid}(\widehat{W}^{2\omega\Delta t})$,

but cells with class rank [math] will undergo four stages:

• a-0.

$\displaystyle\widehat{W}^{\omega\Delta t}=W^{0}+\omega\Delta tR(W^{0})$

• b-0.

$\displaystyle W^{\Delta t}=W^{0}+{\Delta t}R^{Hybrid}(\widehat{W}^{\omega\Delta t})$

• c-0.

$\displaystyle\widehat{W}^{2\omega\Delta t}=W^{\Delta t}+\omega\Delta tR(W^{\Delta t})$

• d-0.

$\displaystyle W^{2\Delta t}=W^{\Delta t}+{\Delta t}R^{Hybrid}(\widehat{W}^{2\omega\Delta t})$.

The different steps for the hybrid time integration with adaptation are quite similar with those presented in Sec. 4.2 for the time update using the full explicit Heun’s scheme.

$\bullet$ Steps 1 and 2:

They areequivalent to those presented for the full explicit scheme in Sec. 4.2

$\bullet$ **Step 3:

**For the cells of class rank [math] sharing a face with a cell of class rank $1$, the update of the solution associated with stage b-0 needs the definition of the flux at the interface $\beta+1/2$,

[TABLE]

$\bullet$ **Step 4:

**It is also equivalent to the one presented for the full explicit scheme in Sec. 4.2. The parabolic interpolation of the residual is performed:

[TABLE]

$\bullet$ **Step 5:

**The approach to time integrate class rank 0 from $\Delta t$ to $2\Delta t$ is quite equivalent to the full explicit time integration. Indeed the flux computed in Eq. (57) is used to compute residual in cell $\beta$ for stages c-0 and d-0. Nevertheless, the following time integration at cell $\beta$ from $\Delta t$ to $2\Delta t$ is chosen :

[TABLE]

This modification with respect to the full explicit Heun’s method is important to ensure space-time conservation of flux from [math] to $2\Delta t$, as shown in Sec. 7.1.

7 Analysis of the conservation and the adaptive time accuracy

7.1 Conservation Property

According to the previous one dimensional configuration presented in Sec. 6.1, the time integration at cell $\beta+1$ of class rank $1$ is formulated as

[TABLE]

As noticed in Sec. 4.3, all negative contributions from the interface $\beta+\frac{1}{2}$ must be recovered in the flux balance of the cell $\beta$.

In order to demonstrate conservation, it is mandatory to integrate the solution in cell $\beta$ until $t=2\Delta t$, with cell $\beta$ belonging to class rank [math]. The time integration leads to:

[TABLE]

in a first phase, and

[TABLE]

in a second phase.

If ${W}^{\Delta t}_{\beta}$ of Eq. (61) is replaced in the second equation of Eq. (62), and reminding that Eq. (57) holds, then

[TABLE]

To conclude, the hybrid time integration in cells $\beta$ and $\beta+1$ from $t=0$ until $t=2\Delta t$ involves the same flux at the interface $\beta+\frac{1}{2}$ and the procedure is conservative.

In the next section, attention is paid on the accuracy of the proposed time-adaptive AION scheme.

7.2 Time Accuracy

The goal is to perform the very same analysis as the one performed for the Heun’s scheme with time adaptation presented in Sec. 5.1. For the conciseness, the full demonstration is not repeated here and only the final form of the error is provided.

Starting again from flux that are $p^{\mbox{\scriptsize th}}$-order accurate in space and from a first-order finite difference approximation of the time derivative, error can be computed for cells of class 0 and 1. Indeed, the numerical error of $W^{2\Delta t}$ at cells of class rank $1$ is:

[TABLE]

and the local error $e({W}^{2\Delta t})$ for cells of class 0 is:

[TABLE]

According to the CFL condition that links $\Delta t$ and $\Delta x$ for a hyperbolic equation ($\Delta x\simeq\Delta t$), the local numerical error is $\mathcal{O}(\Delta t^{2},\Delta x^{p})$ for both class ranks. In the next section, our goal is to recover this theoretical behaviour with numerical simulations.

7.3 Numerical assessment of the theoretical behaviour

The set of equation and the mesh are the same as introduced in Sec. 5.2. In this context, we remind that $y_{k}^{t}$ represents the averaged unknown in the cell $k$ at time $t$. Here, the main change concerns the value of $\omega$ over the computational domain. The parameter $\omega$ is manufactured according to:

[TABLE]

with $\alpha=0.90$, $N$ the number of cells and $j$ the cell index. The hybrid time-adaptive integrator is applied if $0.6<\omega<1$. The value $\omega=0.72$ is chosen arbitrarily for switching between cells of class rank 0 and those of class rank 1 for the time-adaptive AION scheme. For completeness, errors are also provided for Heun’s scheme, time-adaptive Heun’ scheme, AION scheme and time-adaptive AION scheme in Tab. 3. As before, the time adaptive versions of the schemes are referred as +TA. It is numerically demonstrated that the temporal adaptive methods recover the expected orders of accuracy.

The analysis presented in this section gives an information on the local and global accuracy of the scheme, but does not give the two essential behaviours for unsteady simulations: dissipation and dispersion. Section 7.4 is dedicated to the Von Neumann analysis of the time adaptive AION scheme.

7.4 Space-Time von Neumann Analysis

The procedure to perform the von Neumann analysis of the time-adaptive AION scheme follows exactly the same steps as the one for the time-adaptive Heun’s scheme presented in Sec. 5.3. For conciseness, all the details are not provided but the main steps to derive mathematical expressions are summarized.

First, the computational domain introduced in Fig. 12 is considered again, and $\omega$ is set manually according to:

[TABLE]

with $\alpha=0.90$, $N$ the number of cells and $j$ the cell index. The hybrid time-adaptive integrator is applied if $0.6<\omega<1$.

7.4.1 Step UP

In a second step, expressions for the time integration of the different cells are theoretically recovered. Indeed, it leads to the following expressions for cell $\beta$ (step UP) until $t=2\Delta t$:

[TABLE]

Thanks to the hybrid reconstruction for one dimensional advection with a positive advection velocity $c>0$, it comes:

[TABLE]

and the following finite volume formulation in case of hybrid time integration is obtained:

[TABLE]

Here again, the terms $\widehat{y}_{\beta}^{2\omega\Delta t}$ and $(\widehat{\nabla y})^{2\omega\Delta t}_{\beta}$ could be source of instability, as for the time-adaptive Heun’s scheme.

7.4.2 Step DOWN

For the cell $\gamma$ (step DOWN), the hybrid time integration can be formulated as:

[TABLE]

7.4.3 Mathematical expressions

In order to obtain the amplification factor $G$ between $y^{2\Delta t}$ and $y^{0}$, all terms in Eqs. (70) and  (71) must first be expressed according to state $y^{\Delta t}$:

[TABLE]

because of the presence of term $y^{2\Delta t}$ in the right-hand sides of Eqs. (70) and (71). The final expression is finally derived from $y^{\Delta t}$:

[TABLE]

7.4.4 Partial conclusion

The analysis is performed for several values of the parameter $\omega_{j}$ of the hybrid scheme in cell index $j$ since $\omega_{j}$ represents the blending of two schemes. In the following, four values are considered: 0.9, 0.81, 0.72 and 0.65. Figs. 19 and 20 show the dissipation term as a function of $k\Delta x$ for the step DOWN and UP configurations respectively at CFL=0.1. A focus for the low normalised wave numbers ($0<k\Delta x<\pi/6$) shows that amplification occurs for both configurations, something which was not found for the Heun+TA scheme. In addition, amplification becomes greater when the CFL number is increased to 0.3 or 0.6, as shown in Fig. 21 for the step DOWN configuration and in Fig. 22 for the step UP configuration.

Fig. 23 shows that the value of $\omega$ has a minor influence in the dispersion behaviour at CFL=0.1. The same result is obtained at CFL=0.6 in Fig. 24, but the shape of the analysis is more complex, due to the discontinuity in phases. This discontinuity in phase was encountered for the AION scheme in [16] and is a consequence of the specific shape of the transfer function between the parameter $k\,\Delta x$ and the dispersion $\Phi$. Actually, the spectral behaviour for AION+TA and Heun+TA schemes is provided in Figs. 25-27 for two values of the CFL number and for the specific choice $\omega=0.72$. Moreover, the stability limitation of CFL$\leq 0.5$ with Heun time integration is overcome thanks to the hybrid time integration: Fig. 27 shows an amplification of the Heun’s scheme for any wavenumber while amplification only occurs in a short range of wavenumbers for the AION scheme.

7.5 Analysis of $q-$waves

The analysis presented in Sec. 5.4 is now applied to AION+TA configuration, and the effect of hybrid time integration on the occurrence of $q-$waves is explained for configurations UP and DOWN. The analysis follows again the definition of the group velocity and areas of negative group velocity are looked for. The grey zone in Fig. 28 shows the area of negative-$V^{gN}$ group velocity for the AION+TA configuration and steps UP and DOWN. It can be concluded that negative group velocity waves appear essentially for large wavenumbers and the area is larger for the step DOWN configuration than for the step UP and the associated CFL values differ. Indeed, in case of hybrid time integration of Step DOWN, the CFL limit of the negative-$V^{gN}$ group velocity zone is CFL$=0.65$ whereas the CFL limit is $0.3$ in Heun part (Fig. 17). The last question concerns the damping of these $q-$waves using the dissipation of the scheme.

According to Fig. 29, it appears that in case of hybrid time integration, the CFL limit of the negative-$V^{gN}$ zone corresponds to a dissipation $\mu>0.2$ for step UP and $\mu>0.1$ for step DOWN. As a consequence, numerical dissipation can help in attenuating the $q-$waves.

The theoretical analysis provided in this paper is valid for regular grids only and it is of great importance to analyse the time-adaptive AION scheme ability to keep global accuracy for irregular grid. Several configurations are analysed in Sec. 8.

8 Validation

This section is devoted to the validation of the AION+TA scheme using several test cases of increasing complexity, starting from 1D propagation problem to 2D Euler computation.

8.1 Wave propagation problem

The previous theoretical analysis highlighted the fact than Heun+TA and AION+TA do not behave similarly at the interface between temporal classes for regular grid, considering that the theoretical analysis for irregular grid is cumbersome (see Vichnevetsky [26]). Here the goal is to perform the numerical analysis of Heun+TA and AION+TA time integration schemes for an irregular mesh. To do so, a wave packet is propagated in a domain in order to introduce a frequency content inside the computational domain. A non-periodic computational domain of length $L_{x}=270~{}m$ composed of $N=1024$ cells, is initialized with:

[TABLE]

with $K=200$, $f_{e}=1/\pi$ and $x_{c}=90~{}m$. The wave packet is advected at velocity $c=1~{}m/s$. Moreover, while the theoretical analysis was performed using a fixed mesh size, here, two temporal classes of cells with different sizes are introduced. Starting from the space size $\Delta x_{max}$ for the largest cells, the second class with the refined mesh is defined by:

[TABLE]

and $\Delta x_{max}=0.29~{}m$. Step configurations are localised at $x=120.2~{}m$ (step DOWN) and $x=149.5~{}m$ (step UP). Time integration is first performed with the standard Heun+TA scheme until $t=100~{}s$ at the local time step $\Delta t$ in cells of class rank $1$ and $\Delta t/2$ in cells of class rank [math]. Here, $\Delta t$ corresponds to the time step associated with the largest cells, at CFL$=0.6$, which means that:

[TABLE]

Fig. 30 shows, as expected due to the theoretical analysis, that a sinusoidal component of the wave-packet with a certain wavenumber ($p-$waves) is amplified when the wave-packet goes through the step configurations (at $t=25~{}s$ in Fig. 30). A Fast Fourier Transform (FFT) of the numerical solution $y$ is performed to highlight the phenomena, with the sampling frequency equal to $1~{}Hz$.

The main frequency observed from the FFT of $y$ (named $\Psi(y)$ in the following) at time $t=0~{}s$ (Fig. 31) is equal to $0.1~{}Hz$. This frequency is linked to the discretisation of the signal according to $\Delta x_{max}$. The second dominant frequency obtained at $t=25~{}s$ is equal to $0.05~{}Hz$ and corresponds to the discretisation of the signal according to the space step of the finest part of the domain ($\Delta x_{max}/2$). The amplified $p$-waves observed at time $t=25~{}s$ corresponds to the small value of the FFT obtained at the frequency near $0.46~{}Hz$. This frequency is strongly higher than the main frequencies that composed the signal as it is observed in Fig. 30. Hence, Heun+TA is shown to amplify some waves.

Amplification was shown to appear for a computation at CFL=0.6, near the stability limit of Heun’s scheme. If amplification occurs due to the local instability of the space/time scheme, choosing a lower CFL value could reduce the amplified spectrum and in order to confirm our assumption, the very same simulation is now performed at CFL=0.55 until $t=25~{}s$ with the Heun+TA scheme. Fig. 32 focuses on the wave packet advected until $t=25~{}s$. Two computations, with CFL=0.6 and CFL=0.55 are compared and in order to increase readiness, the numerical solution for CFL=0.55 is translated by -0.2 in y-axis. The amplification present at CFL=0.6 disappears at CFL=0.55. The global stability properties of the Heun+TA scheme involve amplification of $p-$waves at CFL=0.6

The same simulation is time-integrated thanks to the AION+TA scheme. For the hybrid part of the AION scheme, the parameter $\omega_{j}$ for cell $j$ is chosen such that:

[TABLE]

with $\alpha=0.90$, and it is also used for the steps UP and DOWN introduced previously.

According to Fig. 33, no amplified $p-$waves appear in the computational domain using the AION+TA time integration. This status is confirmed by the FFT at the frequency $0.46~{}Hz$ (Fig. 34). Hence, even if the AION+TA scheme uses many ingredients of the Heun+TA scheme, it is shown to be more stable and to allow larger CFL values or stable time steps: AION+TA can be seen as an enhancement of Heun+TA time integrator.

The next test case is dedicated to the analysis of the scheme ability to handle compressible effects.

8.2 Sod’s tube

The coupled space/time analysis was always performed using the local regularity of the flow. Here, our goal is to analyse the scheme behaviour on an academic case including discontinuities. The Sod’s tube is an unsteady inviscid one-dimensional test case where the computational domain of length $L_{x}=1~{}m$ is split into two parts separated by a membrane located initially at $x=0.5L_{x}$. The initial flow is defined as:

[TABLE]

where $L$ refers to the left side and $R$ to the right side of the membrane. At $t=0$, the membrane is broken, and waves move inside the computational domain. The solution is composed of a rarefaction wave, a contact discontinuity, and a shock at the final time $t=0.2~{}s$.

The Euler equations are solved using a 1-exact (second-order) space scheme using the Roe approximate Riemann solver. In order to avoid spurious oscillations and to obtain a TVD solution, the minmod slope limiter [20] is used. The computational domain is composed of $N=300$ cells located in regular parts with a uniform mesh size and in irregular parts with non-uniform mesh size such that :

[TABLE]

with $\alpha=0.973$. Time integration is performed using the temporal adaptive procedure, and cells are associated to rank $K$ using:

[TABLE]

with

[TABLE]

Two temporal classes of cells are obtained such as the local time step of the computation is equal to $\Delta t_{max}/2$ in cells of class [math] and $\Delta t_{max}$ in cells of class $1$, where $\Delta t_{max}$ is the maximal allowed time step in the whole domain. For the AION time integration, the parameter $\omega_{j}$ is controlled as:

[TABLE]

For CFL$=0.1$, the Heun+TA scheme is compared to the AION+TA scheme. Figs. 35 and 36 show the density and velocity profiles, with a focus in the region of the rarefaction wave and near the shock. The temporal adaptive approaches (Heun+TA and AION+TA) lead to results in agreement with the theoretical behaviour.

A second set of computations is performed for CFL$=0.45$. Here, the Heun+TA scheme is unstable, and AION+TA results are compared to those of the standard implicit IRK2 scheme (Figs 37 and 38). Both numerical solutions are very close. Paying attention to the velocity and density near the shock, it seems that the AION+TA scheme leads to smoother results than the IRK2 scheme, and overshoots are dissipated.

From these results, it can be concluded that the temporal adaptive approach AION+TA is in agreement with the requirements to handle shock, contact discontinuity and rarefaction waves. Moreover, the AION scheme improves the explicit temporal adaptive method with an enhanced stability.

The next two-dimensional test case is dedicated to the analysis of the AION+TA scheme accuracy.

8.3 Two-dimensional linear advection of an isentropic vortex

The transport of an isentropic vortex, solution of Euler’s equations, is one of the standard test cases of the High Order Workshop since quality of the results is directly linked with scheme accuracy. Indeed, this problem allows to control the capability of the numerical scheme to preserve vorticity in an unsteady simulation and also to estimate the total order of accuracy. The computational domain is composed of a square domain $[-\frac{L}{2},\frac{L}{2}]^{2}$ ($L=0.1$) with periodic boundary conditions. An isentropic vortex defined by its characteristic radius $R$ and strength $\beta$ is imposed on a uniform flow of pressure $P_{\infty}$, temperature $T_{\infty}$ and Mach number $M_{\infty}$. The vortex is initialized in the center of the computational domain $(x_{c},y_{c})=(0,0)$. The initial state is defined by:

[TABLE]

with:

[TABLE]

$R_{gas}=287.15\text{ }J/kg/K$ is the gas constant and a constant ratio of specific heats $\gamma=1.4$ is considered. The isentropic relation leads to the complete set of unknowns:

[TABLE]

The ”fast vortex” test case is considered here and defined by

[TABLE]

The solution is time-marched during three periods inside the periodic box. The computation is performed with a Successive-Correction 2-exact formulation for the spatial scheme (order three) and time integrated by Heun+TA, standard IRK2 and AION+TA schemes. The scheme is not designed to be TVD, but many computations are performed for this case without the need for slope limiters at any order of accuracy, as shown in results of the High Order Workshop. An irregular domain of $260^{2}$ degrees of freedom (DOF) is considered, and the ratio between the largest and the smallest cell sizes is equal to $11$. Reference [16] highlights that the standard Heun’s scheme is not able to perform such test case configuration for a time step $\Delta t$ imposed according to the CFL of the biggest cells contrary to the AION scheme. Here, the computation is designed to obtain four temporal classes of cells, with a $\Delta t_{min}$ of the temporal classes [math] corresponding to CFL$=\frac{U_{\infty}\Delta t_{min}}{\underset{j}{\text{min}}(h_{j})}=0.9$. The parameter $\omega_{j}$ is defined according to cell size as

[TABLE]

with $|\Omega_{j}|$ the volume of cell $j$.

In this case, the proportion of cells with $\omega_{j}<1$ corresponding to hybrid and implicit cells is equal to $73\%$ (see Fig. 39).

According to the pressure and velocity fields in Figs. 40 and 41, it appears that all time integrators have slightly the same properties of dissipation and dispersion. The AION+TA scheme tends to less dissipate than the other second-order time integrators. The temporal adaptive approach keep dispersive and dissipative properties of standard time integrators (Heun and AION) in multi-dimension.

Hence the AION+TA scheme maintain the accuracy of the standard AION scheme.

9 Conclusion

Today industrial complex simulations require at the same time the computation of unsteady effects of turbulence, while maintaining a low CPU cost by optimizing the discretization in the boundary layer and by switching to the RANS model. The work presented in this paper is a first step towards an efficient coupling of RANS and LES simulations in an industrial context. LES and RANS model do not require the same numerical ingredients: LES is generally performed using an explicit time integration in order to easily control the spectral properties of the scheme. In contrary, RANS simulations require to converge fast to the steady solution by means of an unconditionally stable implicit formulation.

In a previous paper, the coupling procedure to handle an explicit time integrator and an implicit formulation in a cell-centered finite volume formulation was presented. The new scheme, called AION, couples the Heun’s scheme and the second-order implicit Runge-Kutta scheme using a transition scheme called hybrid for which a specific definition of the interface flux is required.

Today, standard unsteady industrial simulations with FLUSEPA use a specific time adaptive formulation based on Heun’s scheme. This procedure was not described in details and the current paper focuses first on the description of this scheme currently used in the industrial solver FLUSEPA. This is motivated by the need for a clear description of the scheme and by the need to demonstrate basic mathematical properties. Actually, it is demonstrated that the extended scheme remains second-order accurate, is conservative, but can produce amplification of some waves at the interface between cell classes.

In a second step, the same analysis is extended to the AION scheme coupled with time adaptation. As the time adaptive version of Heun’s scheme, the AION scheme with time adaptation needs to gather cells in temporal classes according to their own maximal (stable) time step and uses sub-cycling for time integration of temporal classes until the biggest time step. Several key ingredients were analysed. First, it was demonstrated that the scheme remains conservative using a space/time analysis. Then, the local time analysis revealed that the AION scheme coupled with time adaptation was second-order accurate on regular grids. Finally, the Fourier analysis of the time adaptive version of AION scheme revealed that dissipation and dispersion behaviours are not really influenced by the time adaptive procedure.

The last part of the paper is devoted to the solutions for a set of canonical applications. One- and two-dimensional simulations with exact solutions confirm numerically a posteriori the theoretical properties demonstrated previously. The procedure is shown capable to handle shock, rarefaction wave and a contact discontinuity.

The next improvement concerns the application of the proposed technique to RANS/LES industrial computations. The smallest cells are time integrated implicitly at high CFL number while the extended temporal adaptive time integration is performed elsewhere. The main corner stone will be to match the hybrid parameter $\omega$ with the shielding function that switches between RANS and LES models.