Thibaut Lunet

Julien Bodart

Serge Gratton

Xavier Vasseur

Time-parallel simulation of the decay of homogeneous turbulence using Parareal with spatial coarsening

Keywords: Direct Numerical Simulation, Explicit time integrator, High-order method in space, Navier-

Direct Numerical Simulation of turbulent flows is a computationally demanding problem that requires efficient parallel algorithms. We investigate the applicability of the time-parallel Parareal algorithm to an instructional case study related to the simulation of the decay of homogeneous isotropic turbulence in three dimensions. We combine a Parareal variant based on explicit time integrators and spatial coarsening with the space-parallel Hybrid Navier-Stokes solver. We analyse the performance of this space-time parallel solver with respect to speedup and quality of the solution. The results are compared with reference data obtained with a classical explicit integration, using an error analysis which relies on the energetic content of the solution. We show that a single Parareal iteration is able to reproduce with high fidelity the main statistical quantities characterizing the turbulent flow field.

Introduction

The numerical solution of time dependent partial differential equations is of interest in many applications in Computational Science and Engineering. The recent advent of computing platforms with more, but not faster, processors currently requires the design of new parallel algorithms that must be able to exploit more concurrency to provide a fast time-to-solution. In this respect, parallel-in-time and space-time methods are considered as promising candidates [START_REF] Dongarra | Applied mathematics research for exascale computing[END_REF]. Indeed, such methods enable the exploitation of possibly substantially more computational resources than purely space-parallel methods with sequential time stepping.

A popular parallel-in-time method based on a decomposition of the space-time domain in time is the Parareal algorithm introduced by Lions, Maday and Turinici [START_REF] Lions | A "Parareal" in time discretization of PDE's[END_REF]. Parareal relies on the availability of a cheap coarse time integrator that provides guesses of the solution at several instants. Given these starting values, a fine time integrator is applied concurrently. The results are then used to propagate a correction to the guesses by using the coarse integrator serially over the time slices. As shown in [START_REF] Gander | Analysis of the Parareal time-parallel time-integration method[END_REF], Parareal can be derived as both a multigrid method in time or as a multiple shooting method along the time axis; see [START_REF] Gander | 50 years of time parallel time integration[END_REF] for further comments and details on the classification of time parallel methods. Due to its non-intrusiveness, Parareal is one of the most widely used time-parallel algorithms. Successful applications have been considered in Computational Fluid Dynamics (CFD) [START_REF] Celledoni | Parallelization in time for thermo-viscoplastic problems in extrusion of aluminium[END_REF][START_REF] Fischer | A Parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], neutron transport [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF], plasma physics [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the Parareal algorithm[END_REF] and skin permeation [START_REF] Kreienbuehl | Numerical simulation of skin transport using Parareal[END_REF] to name a few. While a relative efficiency has been obtained on diffusive problems [START_REF] Fischer | A Parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF][START_REF] Kreienbuehl | Numerical simulation of skin transport using Parareal[END_REF], the applicability of Parareal to hyperbolic or advection-dominated problems is still an open issue; see [START_REF] Ruprecht | Wave propagation characteristics of Parareal[END_REF] for a recent analysis and [START_REF] Bal | Analysis of the convergence and of the stability of the Parareal algorithm for partial differential equations[END_REF][START_REF] Staff | Stability of the Parareal algorithm[END_REF] for early attempts. Although modifications of Parareal for hyperbolic problems have been proposed [START_REF] Dai | Stable Parareal in time method for first-and second-order hyperbolic systems[END_REF][START_REF] Eghbal | Acceleration of unsteady hydrodynamic simulations using the Parareal algorithm[END_REF][START_REF] Farhat | Time-parallel implicit integrators for the near realtime prediction of linear structural dynamics responses[END_REF][START_REF] Gander | Analysis of a Krylov subspace enhanced Parareal algorithm for linear problem[END_REF][START_REF] Ruprecht | Explicit parallel-in-time integration of a linear acoustic-advection system[END_REF], these enhancements generally imply a significant overhead that leads to a degradation of the parallel efficiency. Hence this seems to prevent the use of Parareal for a certain class of important problems in Computational Fluid Dynamics.

In this manuscript we investigate this critical issue by applying Parareal to the direct numerical simulation of complex turbulent compressible flows in three dimensions. We specifically focus on the simulation of the decay of homogeneous isotropic turbulence, a canonical test case known as the simplest configuration which incorporates fundamental and relevant turbulence mechanisms [START_REF] Pope | Turbulent Flows[END_REF]. As identified in [START_REF] Wang | Towards scalable parallel-in-time turbulent flow simulations[END_REF], regarding parallel-in-time integration, interesting features of turbulent flows are their chaotic nature and their strongly unsteady behavior, which typically requires explicit integration methods [START_REF] Choi | Effects of the computational time step on numerical solutions of turbulent flow[END_REF]. This test case incorporates these issues at a very moderate cost, making this model problem an interesting instructional candidate for time-parallel algorithms. To the best of our knowledge, the performance of Parareal on such a test case has not yet been studied. In the following, we concentrate on Parareal with spatial coarsening [START_REF] Fischer | A Parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF] and explicit time integrators. We employ the Hybrid Navier-Stokes solver for the spatial discretization, for which excellent parallel scaling properties have been obtained [START_REF] Bermejo-Moreno | Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion grid points[END_REF]. With this setting, we aim at investigating if this application could benefit from time parallelization. For such a purpose, we perform a numerical study on the role of various parameters on the convergence of Parareal with spatial coarsening. This extensive numerical study is the main contribution of the manuscript.

The manuscript is organized as follows. In Sec. 2 we briefly present the Parareal time-parallel time integration method. In Sec. 3 we describe the set of governing time-dependent partial differential equations and explain the relevance of the canonical test case concerning the direct numerical simulation of turbulent flows. Then we present detailed numerical experiments to understand the convergence of Parareal in Sec. [START_REF] Bermejo-Moreno | Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion grid points[END_REF]. Finally, as a conclusion, we draw first lessons in Sec. 5.

Time parallelization using Parareal

We briefly introduce Parareal [START_REF] Lions | A "Parareal" in time discretization of PDE's[END_REF], a popular method for the time parallel solution of nonlinear initial value problems. Then we describe a variant based on spatial coarsening first proposed in [START_REF] Fischer | A Parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF]. Finally, we include a theoretical model for the expected speedup and parallel efficiency of both algorithms.

General setting

Parareal aims at solving the initial value problem of the form

dU dt = f (U (t), t), U (0) = U 0 , t ∈ [0, T ], (1)
with f :

R p × R + → R p , U (t) ∈ R p , U 0 ∈ R p ,
p being the total number of degrees of freedom and T a positive real value. Here, the problem (1) arises from the spatial discretization of a nonlinear system of partial differential equations (PDEs) through the "method-of-lines" [START_REF] Schiesser | The Numerical Method of Lines: Integration of Partial Differential Equations[END_REF]. We decompose the global time interval

[0, T ] into N uniform time slices [t n-1 , t n ], n = 1, • • • , N
, where N is the number of processes to be considered for the time parallelization only. In the following, we denote by U n the approximation of

U at time t n , i.e. , U n ≈ U (t n ). Let F δt tn-1→tn
(U n-1 ) denote the result of approximately integrating (1) on the time slice [t n-1 , t n ] from a given starting value U n-1 using a fine time integrator F with time increment δ t . Similarly, we introduce a second time integrator G (referred to as the coarse propagator, with time increment ∆ t ), which has to be much cheaper than F in terms of computational time, with possible reduced accuracy. Finally, for ease of exposition, we assume in this analysis that an integer number of both δ t and ∆ t covers a time slice exactly. The prediction step of Parareal consists of computing a first guess of the starting values U 0 n at the beginning of each time slice by

U 0 n = G ∆t tn-1→tn (U 0 n-1 ), U 0 0 = U 0 , (2)
with n = 1, • • • , N . A correction iteration is then applied concurrently on each time slice:

U k n = F δt tn-1→tn (U k-1 n-1 )+ G ∆t tn-1→tn (U k n-1 )-G ∆t tn-1→tn (U k-1 n-1 ), (3)
where U k n denotes the approximation of U at time

t n at the k-th iteration of Parareal (k = 1, • • • , K, n = 1, • • • , N ).
While the application of F can be performed independently for each time slice Parareal remains limited by the sequential nature of the coarse integration (3). Hence, Parareal will bring a reduction of the total computational time with respect to a direct time-serial integration, only if the application of G is cheap enough and if the total number of iterations K of Parareal is small. We recall that Parareal converges after N iterations to the approximation of the exact solution [START_REF] Gander | Nonlinear convergence analysis for the Parareal algorithm[END_REF].

Variant of Parareal based on explicit time integrators and on spatial coarsening

The Parareal parallel in time shooting method [START_REF] Gander | Analysis of the Parareal time-parallel time-integration method[END_REF] is generally used in combination with implicit time integrators; see, e.g. , [START_REF] Croce | Parallel-in-spaceand-time simulation of the three-dimensional, unsteady Navier-Stokes equations for incompressible flow[END_REF][START_REF] Eghbal | Acceleration of unsteady hydrodynamic simulations using the Parareal algorithm[END_REF][START_REF] Kreienbuehl | Numerical simulation of skin transport using Parareal[END_REF][START_REF] Steiner | Convergence of Parareal for the Navier-Stokes equations depending on the Reynolds number[END_REF] for applications related to time-dependent PDEs. In this context procedure, the coarse integrator G is obtained by simply choosing ∆ t > δ t . However, this strategy is usually not applicable when explicit time integrators are favored, since the time step is notably bounded by numerical stability conditions (such as the Courant-Friedrichs-Levy (CFL) condition).

As a cure, Ruprecht and Krause [START_REF] Ruprecht | Explicit parallel-in-time integration of a linear acoustic-advection system[END_REF] have considered a coarse propagator G with lower accuracies in both time and space. Using fewer degrees of freedom at the coarse level is also possible as proposed in [START_REF] Fischer | A Parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF] for the numerical solution of the two-dimensional incompressible Navier-Stokes equations. This has been later investigated in [START_REF] Eghbal | Acceleration of unsteady hydrodynamic simulations using the Parareal algorithm[END_REF] for the Detached Eddy Simulation of the three-dimensional incompressible Navier-Stokes equations and for the simulation of two-dimensional plasma turbulence [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the Parareal algorithm[END_REF], respectively. In this setting, we denote by G the corresponding propagator on the coarse spatial grid (which now involves p degrees of freedom with p < p), while R and I represent the spatial restriction and interpolation operators, respectively. The prediction step of Parareal with spatial coarsening is now

U 0 n = I G ∆t tn-1→tn (R(U 0 n-1 )), U 0 0 = U 0 , (4)
with n = 1, • • • , N . Similarly, the correction iteration of the Parareal algorithm with spatial coarsening then becomes

U k n = F δt tn-1→tn (U k-1 n-1 ) + (5)
I( G ∆t tn-1→tn (R(U k n-1 )) -G ∆t tn-1→tn (R(U k-1 n-1 ))), with k = 1, • • • , K and n = 1, • • • , N .
As noted in [START_REF] Ruprecht | Convergence of Parareal with spatial coarsening[END_REF], the convergence of Parareal with spatial coarsening will not only depend on the fine and coarse time propagators but also on the restriction and interpolation operators. Hence, we expect [START_REF] Berry | Eventbased Parareal: A data-flow based implementation of Parareal[END_REF] to obtain a different convergence rate than [START_REF] Bal | Analysis of the convergence and of the stability of the Parareal algorithm for partial differential equations[END_REF]. In this manuscript, we mostly focus on the variant of Parareal based on (4) and ( 5) with explicit time integration for both the fine and coarse propagators. We next address the expected performance of such a variant.

Expected parallel performance of Parareal

In our setting, parallel-in-time integration is considered as a possibility for additional fine grain parallelism on top of an existing coarse grain spatial decomposition. In a preliminary phase, we have decided to simulate the parallelization in time, whereas the parallelization in space is truly implemented on a distributed memory passing system using the Message Passing Interface (MPI) [START_REF] Gropp | MPI : The complete reference[END_REF]. This allows us to predict at a very moderate cost if the time parallelization can be relevant in our study. Hence, modelling the expected performance of Parareal is of utmost importance. We first analyse the parallel speedup, defined as the ratio of the sequential to the parallel execution time for a given number of processes. As pointed out in [START_REF] Arteaga | A stencil-based implementation of Parareal in the C++ domain specific embedded language STELLA[END_REF][START_REF] Aubanel | Scheduling of tasks in the Parareal algorithm[END_REF][START_REF] Berry | Eventbased Parareal: A data-flow based implementation of Parareal[END_REF][START_REF] Ruprecht | A shared memory implementation of pipelined Parareal[END_REF], the Parareal algorithm is flexible enough to accommodate various implementations based on different programming paradigms. In the modeling, we consider a distributed memory implementation to handle the parallelization in time. We refer the reader to [START_REF] Ruprecht | A shared memory implementation of pipelined Parareal[END_REF] for a discussion and analysis of other strategies. We consider a total of N proc processes for the space-time parallelism with N processes being devoted to the parallelization in time. In this setting, each time slice associated to a spatial subdomain is assigned to a process.

We first consider the standard Parareal algorithm and denote by C F the cost of integrating over a given time slice using F δt tn-1→tn and by C G the corresponding cost when using G ∆t tn-1→tn

. Since explicit timeintegration schemes are used with uniform time steps, we expect both C F and C G to be proportional to the time slice length. Hence C F is equal to T F N δt,F , where T F and N δt,F are the computational time related to the application of the fine integrator over one time step and the number of time steps done for one time slice, respectively. A similar expression can be found for C G (i.e. C G = T G N ∆t,G with similar notation).

As advocated in [START_REF] Ruprecht | A shared memory implementation of pipelined Parareal[END_REF], we concentrate on an efficient implementation of Parareal discussed in both [START_REF] Minion | A hybrid Parareal spectral deferred corrections method[END_REF]Sec. 5] and [START_REF] Aubanel | Scheduling of tasks in the Parareal algorithm[END_REF]Sec. 4] making use of pipelining, i.e. , reducing the costs of the coarse propagation in each correction iteration [START_REF] Bal | Analysis of the convergence and of the stability of the Parareal algorithm for partial differential equations[END_REF] from N C G to C G . The estimate of the theoretical speedup of the pipelined Parareal using K iterations is then given by [START_REF] Minion | A hybrid Parareal spectral deferred corrections method[END_REF] The projected parallel speedup (6) has been derived by neglecting the time spent communicating between each time slice, later referred to as communication in time. For an increased accuracy of the performance model, we have decided to include this cost in the analysis and to propose a modification of pipelined Parareal, which slightly reduces the cost induced by the communications in time. Fig. 1(a) sketches the execution diagram of pipelined Parareal. Since the prediction step (2) is sequential, N -1 communications between time slices are required during this phase (see the thin rectangles on the left of Fig. 1(a)). To avoid this offset, we have considered an implementation shown in Fig. 1(b), in which the solution of the prediction step is computed concurrently on each time slice. This removes unnecessary communications between time slices in the prediction step at the cost of redundancy in the computation. Let C t T denote the cost to communicate a single global (fine) solution from a time slice to the next. Hence, the total cost spent in communications in time for K iterations of Parareal is then given by

S(N ) = N C F N C G + K(C F + C G ) = 1 (1 + K N ) C G C F + K N . ( 6
)
K = 0 K = 1 K = 2 K = 3
C T ,K = C t T K k=1 (N -k) = C t T K(2N -K -1) 2 . (7)
Furthermore, the speedup of our implementation of pipelined Parareal can be obtained as

S T (N ) = 1 (1 + K N ) C G C F + K N + C T ,K N C F . (8)
Let denote by C G the cost of integrating over a given time slice using G ∆t tn-1→tn

. A straightforward adaptation of (8) to the case of pipelined Parareal with spatial coarsening using K iterations yields the following projected speedup

S T (N ) = 1 (1 + K N ) C G + C R + C I C F + K N + C T , K N C F , (9)
where C R (C I ) represents the cost of application of the restriction (interpolation) operator to a vector of appropriate dimension, respectively. Neglecting the communications in time is thus only reasonable if

C T , K N C F (1 + K N ) C G + C R + C I C F + K N , (10)
a condition later discussed in Sec. 4.7.

Finally, we deduce the space-time parallel speedup S S,T (N proc ) as

S S,T (N proc ) = S S (N s ) S T (N ), (11)
where S S (N s ) is the speedup brought by the parallelization in space of the fine solver F on N s processes, i.e. ,

S S (N s ) = T s F ,serial T s F ,Ns . (12)
Hence, we deduce that the space-time parallelization is only viable if

S S,T (N proc ) > S S (N proc ), (13)
a condition later analysed in Sec. 4.2.

We will rely on the expected speedup of the spacetime parallel method S S,T (N proc ) and on the parallel efficiency 1 to predict the parallel performance of our model in Sec. 4.2. We next describe our case study and give details on both the fine and coarse solvers. and the CFD solver We are interested in the simulation of complex turbulent flows and focus on a canonical but relevant test case, well known from the CFD community. The decay of homogeneous isotropic turbulence (HIT) has been studied by many authors from moderate [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF] to very large scale problems [START_REF] Kaneda | Energy dissipation rate and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box[END_REF]. This case may be seen as the simplest configuration which incorporates fundamental and relevant turbulence mechanisms [START_REF] Pope | Turbulent Flows[END_REF]. In the context of Direct Numerical Simulation (DNS), the problem size is essentially driven by the range of length scales to resolve, which sets the required grid resolution. Specifically, increasing the Reynolds number leads to a wider spectrum of length scales to resolve, but also improves the relevance of the simulation regarding high Reynolds, industrial configurations. In this setup, synthetically generated turbulence decays under the action of dissipation mechanisms. A description of the flow field during the simulation is given in Fig. 2 in terms of turbulent kinetic energy, which highlights the evolution of energy-carrying turbulent eddies and the strong nonlinearity of turbulent flows.

Regarding parallel in time integration algorithms, interesting features of turbulent flows are i) their chaotic nature ii) their strongly unsteady behavior, which typically requires explicit integration [START_REF] Choi | Effects of the computational time step on numerical solutions of turbulent flow[END_REF]. Interestingly, the decay of HIT incorporates these issues at a moderate cost, making this model problem a challenging candidate for time-parallel algorithms. 2 ).

Governing equations

The compressible Navier-Stokes equations read:

∂ρ ∂t +∇.(ρu) = 0, (14)
∂ρu ∂t +∇.(ρuu + pδ) = ∇.τ , (15) ∂ρe
T ∂t +∇.(u(ρe T + p)) = ∇.(u.τ -q), ( 16
)
where ρ is the density, u = (u, v, w) the velocity vector, p the pressure, δ the unity tensor, e T = e + u.u 2 the total energy, e the specific internal energy, τ the viscous stress tensor and q the heat flux. e, p and ρ are linked by the equation of state, obtained for a calorically perfect gas:

p = (γ -1)ρe, (17)
with γ the specific heat coefficient ratio, whereas τ and q are modeled via classical Stokes' hypothesis:

τ = 2µS - 2 3 µ(∇.u)δ, (18)
q = -k c ∇T, (19)
with S = ∇u + (∇u) T /2 the strain-rate tensor, µ the dynamic viscosity, k c the thermal conductivity, the temperature T = (γ -1)e/R and R the gas constant. Temperature dependency for µ is accounted for through a power-law assumption:

µ/µ ref = (T /T ref ) 3/4 , (20)
and the ratio between µ and k c is set constant through a constant Prandtl number, which closes the set of equations.

Decay of homogeneous isotropic turbulence

The decay of homogeneous isotropic turbulence is studied using a uniform grid on a 3D periodical box of size L. Note that all the flow quantities discussed in this section are time dependent unless otherwise specified.

In particular the 0 subscript refers to the initial state at t = 0. Because the problem is homogeneous, the statistical evolution of HIT decay can be reduced to a temporal evolution using the following spatial averaging:

f = 1 L 3 V f (x, y, z) dxdydz. ( 21
)
Initial conditions may be obtained by setting the flow fields (ρ 0 , u 0 , T 0 ), and parameters R, γ, P r, µ ref and T ref . These fields are constructed through a random process that builds a synthetic turbulent flow field following a prescribed energy spectrum, mostly concentrated at large scales. The detailed methodology associated with the construction of the initial turbulent field is described in Sec. 3.5 and Ap. A of [START_REF] Johnsen | Assessment of highresolution methods for numerical simulations of compressible turbulence with shock waves[END_REF], respectively. This procedure sets the flow conditions, essentially characterized by two non-dimensional, time evolving quantities which drive the turbulence evolution: the Taylor scale Reynolds number2 ,

Re λ = ρ u λ µ , (22)
and the turbulent Mach number:

M t = √ 3 u c , (23)
where u = ( uu + vv + ww )/3, and c = √ γRT is the speed of sound. λ is the (transverse) Taylor micro-scale, and provides a length scale associated with intermediary-sized eddies of the flow field 3 . Re λ sets the width of the turbulent kinetic energy spectrum, and M t indicates the influence of compressibility effects, that may trigger discontinuities (e.g. shocklets) in the solution for high M t . The turbulent kinetic energy k e and the dissipation , are defined as follows:

k e = 3 2 u 2 , = µ ρ S.S . (24)
Because the footprint of k e and is mostly located at large and small scales, respectively, these two quantities are interesting indicators to track the turbulence behavior during the HIT simulation. We can derive a transport equation for k e which reduces, using the 3D homogeneous hypothesis and thus isotropy, to:

dk e dt = -. (25)
It has been observed [START_REF] Pope | Turbulent Flows[END_REF] that, after a short transient phase, both values exhibit a power-law decay, as shown in Fig. 3. Physically, it illustrates the energy cascade of the turbulence: energy from the large scale eddies is continuously transferred to the small scale eddies until the latter dissipate through molecular viscosity. The Kolmogorov length scale characterizes the smallest, dissipative scales:

η = µ 3 ρ 3 1/4 , (26)
which sets the smallest grid cell size to achieve direct numerical simulation. Smallest length scales must be accurately resolved on the spatial grid. Considering the largest wavenumber that can be represented on a given spatial uniform Cartesian mesh of size N 3 L (see [START_REF] Pope | Turbulent Flows[END_REF]Chap. 9]):

κ max (N L ) = πN L L , (27)
we can derive a condition consistent with the correct resolution of eddies of size η:

∀t, η κ max (N L ) ≥ α, (28)
where α is a coefficient dependent of the solver and specifically the space discretization scheme.

Massively parallel Navier Stokes solver Hybrid

We use the compressible structured solver Hybrid, developed in C++ by [START_REF] Larsson | Effect of numerical dissipation on the predicted spectra for compressible turbulence[END_REF], which aims at studying fundamental turbulence problems, such as shock-turbulence interaction [START_REF] Larsson | Reynoldsand Mach-number effects in canonical shock-turbulence interaction[END_REF]. In the absence of discontinuities, the code uses a centered 6 th -order finite difference scheme, while time-integration is performed thanks to a 4 th -order explicit Runge-Kutta method (RK4). Full details on space and time discretizations are given in [START_REF] Johnsen | Assessment of highresolution methods for numerical simulations of compressible turbulence with shock waves[END_REF]. The Hybrid solver uses MPI-based parallelism which relies on space decomposition. It shows very good weak scaling results when using up to 2 million cores [START_REF] Bermejo-Moreno | Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion grid points[END_REF]. The adopted numerical methodology developed in Hybrid (structured mesh, explicit space-time discretization) also exhibits excellent strong scaling properties, if the number of cells N min cpc of the spatial subdomain associated with one MPI process (and thus a single CPU core) exceeds a few hundreds. The exact value of N min cpc to maintain parallel efficiency is necessarily architecture dependent, and sets the lower bound of the operating limit for a sole and efficient space parallelization. Given its moderate efficiency, time parallelization should be used to extend the number of usable core for efficient computations, thus keeping N cpc > N min cpc . Once the simulation is properly initialized, an appropriate choice of the time step needs to be considered to minimize the computational cost of a direct numerical simulation. Although implicit time stepping is usually not rewarding because of the limited accuracy observed at large time steps [START_REF] Choi | Effects of the computational time step on numerical solutions of turbulent flow[END_REF], smallest physical time scales are still significantly larger than the maximum time step δt CF L authorized by stability considerations. Hence, cost efficient simulations must be carried out with δt ≈ δt CF L . For this specific solver, we have estimated the CF L number as 1.79. This value indeed corresponds to the limit of linear stability of the 6 th -order centered finite difference scheme associated with the RK4 time-integration method.

Transfer operators for Parareal with spatial coarsening

We detail next the transfer operators used in the variant of Parareal based on spatial coarsening as briefly introduced in Section 2.2. Since the space discretization in Hybrid is based on high-order finite difference schemes on a Cartesian structured mesh, geometric coarsening in space (known as vertex-centered coarsening in multigrid [START_REF] Trottenberg | Multigrid[END_REF]) is adopted to construct the coarse level of the two-grid hierarchy. In this setting, the fine mesh is coarsened by simply considering every other node in the three spatial directions, thus leading to a cell size ratio of 2 in each direction between the coarse and the fine mesh.

As restriction operator R, we have considered both the injection and the three-dimensional full-weighting operator [43, Sec. 2.9]. Due to this coarsening in space, the time step on the coarse grid is simply set as

∆ t = 2δ t ,
to keep the same CFL condition on both fine and coarse grids. Concerning I, we have selected various interpolations (trilinear [43, Sec. 2.9], tricubic [START_REF] Lekien | Tricubic interpolation in three dimensions[END_REF], 7 th -order or Fourier transformation based [START_REF] Russell | Optimised three-dimensional Fourier interpolation: An analysis of techniques and application to a linear-scaling density functional theory code[END_REF]). The 7 th -order interpolation is based on tensor products of one-dimensional operator, which makes it simultaneously computationally cheap and implementation friendly. The influence of both restriction and interpolation operators on the convergence of Parareal with spatial coarsening as well as their costs will be numerically investigated in Sec.s 4.4 and 4.7, respectively.

Direct Numerical Simulation methodology

In the following subsection, we develop the methodology to perform the HIT Direct Numerical Simulation, with the objective to propose a systematic and reproducible framework for these simulations. The Navier-Stokes equations ( 14), ( 15) and ( 16) are solved following a direct numerical simulation framework. Thus, the grid requirements need to be adjusted with the minimum length scale η, which sets the ratio between the largest scales (i.e. the box size L) and the smallest scales. Recent DNS of HIT with spectral space discretization used α = 1 [START_REF] Kaneda | Energy dissipation rate and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box[END_REF]. Taking into account that Hybrid uses high-order centered space discretization methods, we choose a more restrictive α = 1.5 criterion in [START_REF] Lions | A "Parareal" in time discretization of PDE's[END_REF] for the simulation with F, consistent with modified wave numbers expression for 6 th -order schemes in space (e.g. [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] ). This resolution criterion needs to be satisfied at all times during the decay of turbulence. As the flow field undergoes a transient during which it evolves from synthetic to physical turbulence, with a characteristic "eddy turnover" time scale k e / , we specifically enforce this criterion after a transient t φ , i.e. using η = η(t φ ) in [START_REF] Lions | A "Parareal" in time discretization of PDE's[END_REF].

In order to propose an estimation of t φ , we first extract the spectral density of energy E(κ) of the turbulent flow field. It represents the energy content associated with the norm κ of the related wave vectors, and is interpreted as so called turbulent eddies of corresponding length scale in the study of turbulent flows. Hence,

k e = κ E(κ). ( 29
)
Details on the computation of E(κ) for HIT are provided in Ap. A. A representation of energy spectra is given in Fig. 4 at various snapshots of the HIT decay. The spectrum of the synthetic initial solution has an energy peak located at the large scales, consistent with the initialization process. After some time, the turbulent energy cascade takes place, which results in k e being spread over every length scales of the flow field.

We estimate the end of this transient t φ as the time for which the spectral density of energy starts decaying at all scales:

t φ = min {t > 0 | ∀κ, δE(κ) < 0} , (30)
where δE(κ) is the growth rate of E(κ). This criterion is useful in finding a systematic condition to determine the maximum Re λ0 that can be simulated using a mesh of given size N 3 L . We choose to define it as the maximum value of Re λ0 that satisfies condition (28) at t φ , that is:

Re λ0,max = max {Re λ0 | ∀t > t φ , η κ max (N L ) ≥ α} . (31)
A trial and error procedure is therefore required to obtain the desired resolved state at time t φ and results in a well identified initial state, fully characterized by Re λ0 .

To complete the methodology, an error metric is defined to estimate the accuracy of the parallel time integration. As a basis in the application of Parareal to turbulent flows we analyse the physical relevance of the solution rather than the strict convergence of the algorithm. Because the complex flow dynamics of a turbulent flow is characterized in a statistical sense by a transfer from the large to the small scales, the quality of the solution is evaluated through the spectral density of energy. It may reveal undesired effects such as nonphysical dissipation or dispersion promoted by the numerical scheme. The accuracy of the solution is measured using an error metric developed in Section 4.3, which relies on the energy content of the flow field and characterizes the ability of Parareal to capture the flow dynamics.

Numerical results

This section provides a detailed numerical study of the different factors influencing the parallel performance of Parareal with spatial coarsening, when considering our test case in Computational Fluid Dynamics. In Sec. 4.2, we compare the pure spatial and space-time parallelizations in terms of speedup and efficiency. This allows us to determine situations for which space-time parallelization is meaningful. Then, in this setting, we examine the numerical quality of the solutions provided by Parareal in Sec. 4.3. In Sec. 4.4, 4.5 and 4.6 we investigate the influence of several parameters on the quality of the solution. Finally, in Sec. 4.7, we continue investigating the parallel efficiency of Parareal with spatial coarsening, by including in the analysis the costs of the transfer operators and of the time communications.

Methodology

The numerical simulations were performed on the EOS supercomputer at CALMIP, Toulouse, France. This platform is equipped with 612 compute nodes, each node hosting two 10 core Intel Ivy Bridge chips (that run at 2.8 GHz), and 64 GB of system memory. Each socket was equipped with 25 MB of cache memory. EOS nodes are connected by a nonblocking fat tree network, with a network bandwidth of 6.89 Gbytes/s for sending and receiving data. The code was compiled using the Intel 14.0.2.144 compiler, with Intel MPI 4.1.3.049 library.

As pointed out in Sec. 2.3, the parallel performance of Parareal is simulated. More precisely, the computations related to all time slices in Parareal are performed sequentially in time, whereas the applications of both F and G through the Hybrid code are parallel in space. To estimate the cost to communicate a single global solution C t T in Parareal, we have adopted the following simple strategy. First, we recreate exactly the same spatial decomposition on N proc /N processes as in Hybrid. Secondly, we perform standard MPI communications of all the relevant physical fields from one time slice to the next over the N processes and measure the corresponding wall-clock times on the EOS supercomputer. This procedure is then repeated 20 times to provide a meaningful estimation, C t T being then obtained as the maximum wall-clock time. We refer the reader to the C++ source code 4 for further details. A similar methodology is applied to obtain estimates of C F and C Ĝ .

In the following sections, unless stated otherwise, we consider the direct numerical simulation of the decay of homogenous isotropic turbulence at Re λ = 46 on a Cartesian mesh of size 80 3 (i.e. N L = 80). We have estimated the initial CFL number as 1.79 for each simulation. Due to the energy decay, the CFL number is found to slightly decrease as time increases.

Preliminary analysis of the Parareal parallel performance

We address the question of the parallel performance of Parareal based on spatial coarsening. An important point is to determine whether space-time parallelization can be more appropriate than plain spatial parallelization in our application. To this end, we operate the solver with N cpc < N min cpc to deliberately reach the limit of the spatial decomposition. Given F and G introduced in Section 3.3, we consider the influence of both N and K (number of time slices and number of iterations, respectively) on the speedup of Parareal with spatial coarsening. For ease of exposition, we first neglect the costs related to both the transfer operators (C R and C I ) and the communications in time (C T, K ) in [START_REF] Dai | Stable Parareal in time method for first-and second-order hyperbolic systems[END_REF]. A refined analysis will be proposed later in Sec. 4.7.

The strong scalability of F and G is shown in Fig. 5 (top). The fine solver F exhibits good scalability properties for a number of processes up to 160. Nevertheless, for a larger number of processes, the performance starts to saturate. Not surprisingly, we also observe a rather quick deterioration of the parallel performance of the coarse solver G. It is important to remark at this stage that this earlier decrease in performance of the coarse solver will bound the cost ratio C F /C G at a lower value than the expected one 5 (which is equal to 16, assuming a perfect strong scaling and excluding communication costs).

Given a number of processes in space fixed to N proc /N = 160, we consider the additional parallelization in time with 2, 4, 8 and 16 time slices, respectively. The simulated speedups and parallel efficiencies obtained after one or two iterations of Parareal with spatial coarsening are then provided in the 4 available at https://gitlab.com/tlunet/parallel-in-time 5 This value of the theoretical ratio of coarse-to-fine execution is explained by the coarsening in space by a factor of 2 in each direction and by the choice of the coarse time step detailed in Sec. 3.4. Coarse solver, N = 1 to the pure spatial parallelization, whatever the number of time slices. This is a rather satisfactory result. However, for two or more iterations of Parareal with spatial coarsening, the success is limited. Indeed, for N proc = 320 or N proc = 640, the speedup of the spacetime parallel algorithm is found to be lower than the speedup obtained with exclusive space parallelization. Not surprisingly, the parallel efficiency is maximized if the number of Parareal iterations is limited. Interestingly, the gap of parallel efficiency between K = 1 and K = 2 reduces as the number of time slices N is increased. This trend is explained because the reduced efficiency of the coarse solver is balanced by a larger number of fine solver computations as K increases.

K = 1, C F /CĜ = 16.0 K = 1, C F /CĜ = 9.4 K = 2, C F /CĜ = 16.0 K = 2, C F /CĜ = 9.4 K = 1, C F /CĜ = ∞ K = 2, C F /CĜ = ∞
However, because the considered efficiency levels are rather low, we mainly favor the application of a single iteration of Parareal with spatial coarsening in the rest of this study 6 . We now investigate the influence of the N proc /N parameter on the performance of the combined spacetime parallelization in our application. Hence, we introduce the relative gain of the space-time parallelization as

σ S,T = S S,T (N proc ) -S S (N proc ) S S (N proc ) . ( 32
)
A positive value of σ S,T thus indicates that the spacetime parallelization is worth considering. We denote by σ S,T the relative gain of the space-time parallelization when neglecting the costs of communication in time and transfer operators. Although σ S,T is an ideal value, it is worth noting that σ S,T ≈ σ S,T for a sufficiently long time slice, i.e. if C T ,K , C R and C I are negligible compared to C F and C Ĝ .

Table 1: Relative gain of the space-time parallelization σ S,T after one iteration of Parareal with spatial coarsening (neglecting the costs of communication in time and transfer operators). N proc denotes the total number of processes with N proc /N processes being assigned to the parallelization in space only. The numbers of time slices considered here are 2, 4, 8 and 16, respectively.

N proc 160 320 640 1280

N proc /N = 80 -18% -3.9% +10% +49% N proc /N = 160 - +0.82% +24% +82% N proc /N = 320 - - -1.1% +54%
Tab. 1 collects the values of σ S,T versus N proc /N , when a single iteration of Parareal with spatial coarsening is performed. For a fixed number of processes devoted to the spatial parallelization, increasing the number of time slices (N ) does improve σ S,T as expected. This behavior is indeed in agreement with Fig. 5 (top). Tab. 1 also reveals that the parallelization in space and in time is worth considering when the total number of processes is large. In this situation, we rather favor the case of a low number of time slices (see the bold values in Tab. 1), since the convergence of Parareal is reached in at most N iterations. Hence, in what follows, we choose N = 4 and study the quality of the solution of the first iterations of Parareal on a large number of processors (N proc = 640 or N proc = 1280). This is next investigated in Sec. 4.3.

Energy spectra after successive iterations of Parareal with spatial coarsening

We here analyze the flow solution obtained after successive iterations of Parareal with spatial coarsening. In particular, we focus on the energy spectrum of the resulting turbulent flow. As described above, we use trilinear interpolation and injection as transfer operators, and N ∆t, Ĝ = 20. The simulation of the HIT decay is solved between physical times t φ = 1.4T and t end = 4.4T , which corresponds to the relevant part of the energy decay (see Fig. 3). The energy spectrum is represented on Fig. 6 for Ĝ (coarse integrator without interpolation), F (fine integrator, reference solution) and the first iterations of Parareal with spatial coarsening, respectively. Because of the insufficient spatial resolution for Ĝ, an over-estimation of the energy on the small scales around κ = 0.5κ max . This is a classical observation of unresolved direct simulation that results from a partial cut off of the physical dissipation mechanisms operating at the smallest scales. The effects of the interpolation are directly highlighted in the case K = 0 for which the middle range scales are significantly diffused, while a peak of energy is observed for the smallest scales. The first iteration ( K = 1) results in a dampening of this high frequency peak, along with a significant improvement of the middle range scales. Interestingly the second iteration ( K = 2) does not bring significant improvements, reinforcing the idea of using a single iteration of Parareal with spatial coarsening.

K = 0 K = 1 K = 2
To better quantify the quality of the solution, we extract the relative error of the energy spectrum e rel , defined as:

e rel (κ) = Ẽ(κ) -E F (κ) E F (κ) (33)
where Ẽ(κ) corresponds to the spectrum of an approximate solution, ( Ĝ or Parareal iteration) and E F (κ) the spectrum of the reference solution computed with F solely. While a negative value for e rel (κ) indicates energy damping, a positive value indicates energy amplification. We add quantification on the spectrum error by looking at the relative error on quantities integrated over all length scales, the energy:

e ke = ke -k e,F k e,F , (34)
and the dissipation:

e = ˜ -F F . (35)
These two scalars are representative of statistical errors at two levels: e ke characterizes the large length scales, while e is more a measure of the behaviour at small length scales. The three error indicators e rel , e ke and e will be combined to investigate the influence of several parameters on the quality of the solution: the transfer operators in Sec. 4.4, the Reynolds number in Sec. 4.5 and the time slice length in Sec. 4.6.

Influence of the transfer operators

The influence of the restriction operator on the relative energy error is shown in Fig. 7. In agreement with the developments above, we consider the solution after one iteration of Parareal to provide a possible comparison. Whatever the restriction operators, the absolute value of the relative energy error is found to grow in the small scale regime. Moreover, the curves related to the injection and the full-weighting operators look very similar. Hence, we rather consider the injection as restriction operator in the rest of the manuscript, to minimize C R . The influence of the interpolation operator on the relative energy error is shown in Fig. 8. The corresponding relative error indicators on the turbulent kinetic energy and on the dissipation respectively are given in Tab. 2. Fig. 8 reveals that high-order interpolation operators must rather favored to minimize the relative error in energy. This is also confirmed in Tab. 2. Hence, in the following, we consider the 7 thorder interpolation operator as a standard choice.

Influence of the Reynolds number

We consider the simulation of the decay of homogeneous isotropic turbulence on two finer grids (160 3 and 320 3 , respectively) associated with larger Reynolds numbers to analyse the influence of the Reynolds number on the quality of the solution. From a physical standpoint, the major effect on the simulation is a widening of the range of length scales carried by the turbulent flow.

The influence of the Reynolds number Re λ on the relative error on the energy is shown in Fig. 9. The three curves exhibit a very similar behavior at all scales for each considered Reynolds number. The corresponding relative errors on the turbulent kinetic en- We thus investigate the influence of N ∆t, Ĝ on the energy relative error after one iteration of Parareal with spatial coarsening. Since changing this parameter will also modify the physical time t end , we have decided to restart the time parallel algorithm after its first iteration to complete the full simulation with t end matching the other choices of N ∆t, Ĝ . We consider a simulation at Reynolds number Re λ0 = 322 with several values of N ∆t, Ĝ ∈ {10, 20, 40}, which implies restarting Parareal with spatial coarsening 3, 1 and 0 times, respectively. Hybrid was ran using the same CFL constraint and a number of time slices of N = 4, leading to t end = 4.8T .

The main effects of the variation of N ∆t, Ĝ on e rel are represented in Fig. 10. We observe a reduction of the dampening error at the small scales. We attribute this effect to the dissipation of the small scales initially affected by injection/interpolation steps. Increasing N ∆t, Ĝ does not significantly deteriorate e , which remains almost independent of N ∆t, Ĝ , although we observe a slight dampening appearing for the middle range, more energetic length scales (κ ≈ 0.4κ max ). On the other hand, the latter amplification has an expected detrimental effect on the level of e ke which increases with N ∆t, Ĝ . Results are summarized in Tab. 4, together with the results obtained with the coarse integrator operated solely as a reference. Note that the apparently small error observed with the coarse inte- grator benefits from an error compensation with the adopted metric.

Refined analysis of the speedup of Parareal

Finally, based on (9), we include the costs of both the communications in time (C t τ ) and the transfer operators (C I and C R ) in the analysis of the speedup. We consider the configuration at Re λ0 = 46 (N L = 80) with the 7 th -order interpolation and the injection as transfer operators.

The cost related to the restriction operator can be neglected since no floating operations are involved, i.e. , C R = 0. Furthermore, we can relate the cost of the interpolation to the cost of performing one single time-step of the coarse solver. Indeed the onedimensional formulas for either the 7 th -order interpolation or the 6 th -order finite difference discretization in Hybrid involve a similar complexity (15 and 11 floating point operations are required, respectively).

A more precise estimation of the complete computational complexity 7 then leads to

C I ≤ 15 39 × 15 11 C Ĝ N ∆t, Ĝ . (36)
In the following, we consider the upper bound in [START_REF] Ruprecht | Wave propagation characteristics of Parareal[END_REF] as an accurate estimation of C I . The relative gain of the space-time parallelization σ S,T (see [START_REF] Reynolds-Barredo | Mechanisms for the convergence of time-parallelized, Parareal turbulent plasma simulations[END_REF]) is then given in Tab. 5. Tab. 5 reveals that the effective relative gains σ S,T are relatively close to the ideal values σ S,T discussed in Sec. 4.2. The effective gain can be indeed significant (especially when considering a large number of processors), which is a satisfactory behavior. Finally, comparing cases (a) and (b) in Tab. 5 reveals that the influence of C I is found to be marginal.

Conclusions

This manuscript investigates the applicability of Parareal with spatial coarsening to a relevant test case in Computational Fluid Dynamics related to the simulation of the decay of homogeneous isotropic turbulence. The time parallel simulation of such flows does actually lead to an instructional test case, since their main interesting features are both their chaotic nature and their strongly unsteady behavior. In a first phase, we have decided to simulate the parallelization in time, whereas the parallelization in space is truly implemented on a distributed memory passing system. Explicit Runge-Kutta time integration methods 7 In Hybrid the computation of the right-hand side of the Navier-Stokes equations involves 39 gradient evaluations (3 for (14), 18 for (15), 6 for (16), 9 for (18) and 3 for [START_REF] Gropp | MPI : The complete reference[END_REF], respectively). Since the RK4 time integration method requires 4 evaluations of right-hand sides, the total amount of gradient evaluations for one time step of the coarse solver is then 156. Finally, 60 stencil evaluations per coarse point are required for the interpolation of the 5 different fields. and high-order finite difference schemes are used for the temporal and spatial discretizations of the Navier-Stokes equations.

A methodology related to the computation of the energy spectrum has been proposed to assess the numerical quality of the iterative solution provided by the time parallel algorithm. Based on this analysis, we have found that the solution after a single iteration of Parareal with spatial coarsening was physically relevant, provided that a high-order interpolation operator in space is employed. In this setting, the extensive numerical experiments clearly illustrate the possible benefits of using parallelization in time. This rather encouraging result from a physical point of view needs of course to be confirmed by a detailed convergence study. This is indeed an important research direction that we are currently considering.

We are deeply convinced that this test case can serve as a relevant benchmark for time parallel methods. Hence, to propose a reproducible framework, we have carefully described the complete methodology and plan to make the simulation code freely available to the community, if a significative interest appears. Indeed, an instructive next step would be then to assess the performance of current or emerging time parallel algorithms on this configuration.

A Computation of the energy spectrum

We briefly describe the computation of the energy spectrum. More details can be found in the literature (see, e.g. , [20, Ap. B]). We denote by û(κ x , κ y , κ z ) the Discrete Fourier Transform (DFT) of u(x, y, z), using 1 (1/N 3 L ) to normalize the direct (inverse) DFT, respectively. We define κ x , κ y and κ z the discrete one-dimensional wavenumbers as:

κ [.] = n 2π L with n ∈ Z | - N L 2 + 1 ≤ n ≤ N L 2 , (37)
and the associated vector norm:

κ = κ = κ 2 x + κ 2 y + κ 2 z . (38)
The 3D Fourier space is discretized into three-dimensional shells of thickness dκ (usually 2πdκ L = 1). To deduce the energy corresponding to the wavenumber κ i , a bin count is then performed using an integer number of shells: (39)

E(κ i ) = 1 2N 6 κ∈[κ i -dκ 2 ,κ i + dκ

Fig. 1 :

 1 Fig. 1: Execution diagram of Parareal with pipelining (top) and with improved pipelining (bottom). Three iterations of Parareal (K = 3) on four time slices (N = 4) are considered. The end of each Parareal iteration is indicated by a dotted line. Time-slices are numbered from 0 to 3, each line corresponds to one process of Parareal.

Fig. 2 :

 2 Fig. 2: Plane cuts of local turbulent kinetic energy k e,loc /k e,0 for Re λ,0 = 322 (k e,loc = u.u 2 ).

Fig. 3 :

 3 Fig. 3: Time-evolution of k e /k e,0 and /k e,0 , with T = (λ/u ) t=0 , N L = 80, Re λ 0 = 46.

Fig. 4 :

 4 Fig. 4: Energy spectra during the decay of HIT, with T = (λ/u ) t=0 , N L = 80, Re λ 0 = 46.

SpeedupFine solver, N = 1 1 K = 1 , 4 20Efficiency

 1114 Coarse solver, N = C F /CĜ = 16.0 K = 1, C F /CĜ = 9.4 K = 2, C F /CĜ = 16.0 K = 2, C F /CĜ = 9.Fine solver, N = 1

Fig. 5 :

 5 Fig. 5: Simulated speedup (top) and parallel efficiency (bottom) of Hybrid (fine and coarse solvers), and of Parareal with spatial coarsening after ( K = 1, 2) iterations.

Fig. 5 .

 5 Fig.5. The inherent limits of Parareal are indicated in Fig.5b. They are obtained by considering an infinite cost ratio C F /C G . Performing a single iteration of Parareal with spatial coarsening does lead to an increased speedup and parallel efficiency with respect to the pure spatial parallelization, whatever the number of time slices. This is a rather satisfactory result. However, for two or more iterations of Parareal with spatial coarsening, the success is limited. Indeed, for N proc = 320 or N proc = 640, the speedup of the spacetime parallel algorithm is found to be lower than the speedup obtained with exclusive space parallelization. Not surprisingly, the parallel efficiency is maximized if the number of Parareal iterations is limited. Interestingly, the gap of parallel efficiency between K = 1 and K = 2 reduces as the number of time slices N is increased. This trend is explained because the reduced efficiency of the coarse solver is balanced by a larger number of fine solver computations as K increases.

Fig. 6 :

 6 Fig. 6: Energy spectra after K iterations of Parareal with spatial coarsening at t end (N = 4, N ∆ t , Ĝ = 20, trilinear interpolation and injection as transfer operators).

7 :Fig. 8 :

 78 Fig. 8: Influence of the interpolation operator on the relative energy error after one iteration of Parareal with spatial coarsening (N = 4, N ∆ t , Ĝ = 20, injection as restriction operator).

Fig. 9 :

 9 Fig. 9: Influence of the Reynolds number Re λ on the relative energy error e rel after one iteration of Parareal with spatial coarsening (N = 4, N ∆ t , Ĝ = 20, injection and 7 th -order interpolation as transfer operators).

Fig. 10 :

 10 Fig.10: Influence of N ∆ t , Ĝ on the relative error on the energy spectrum after one iteration of Parareal with spatial coarsening (N = 4, injection and 7 th -order interpolation as transfer operators).

2 [|

 2 û(κ x , κ y , κ z )| 2 .

Table 2 :

 2 Influence of the interpolation operator on the relative errors on k e and after one iteration of Parareal with spatial coarsening (N = 4, N ∆ t , Ĝ = 20, injection as restriction operator).

		Interpolation	e k e		e
		Trilinear	-6.6%	-17%
		Tricubic	-1.5%	-5.5%
		7 th -order -0.85% -2.8%
		Fourier -0.79% -2.1%
		Ĝ alone -0.83% -2.7%
	0.0	0.2	0.4	/ max	0.6	0.8	1.0

Table 3 :

 3 Influence of the Reynolds number Re λ on the relative error on k e and after one iteration of Parareal with spatial coarsening (N = 4, N ∆ t , Ĝ = 20, injection and 7 thorder interpolation as transfer operators). Re λ φ (Re λ E ) is the Reynolds number at the beginning (end) of the parallel in time simulation.N L Re λ 0 t φ /T t end /T Re λ φ Re λ E Parareal to handle correctly the dissipative scales of the turbulent flow. The grid size being designed to specifically handle smallest length scales, this constraint does not evolve significantly relatively to the cell size, hence with an increase of the Reynolds number. These encouraging results will be the subject of a forthcoming study. It suggests similar behavior of Parareal in large scale simulations characterized with significantly higher Reynolds numbers (e.g.[START_REF] Kaneda | Energy dissipation rate and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box[END_REF]).4.6 Influence of the time slice lengthThe increase of the Reynolds number previously developed is also useful to test the influence of the time slice size in the simulation, as investigated in[START_REF] Eghbal | Acceleration of unsteady hydrodynamic simulations using the Parareal algorithm[END_REF] Sec.4], [14, Sec 5.2] and[START_REF] Reynolds-Barredo | Mechanisms for the convergence of time-parallelized, Parareal turbulent plasma simulations[END_REF] Sec 4.5]. This is indeed a critical parameter regarding the parallel efficiency of the algorithm as mentioned in Section 2.3. Indeed, for small Reynolds numbers, the required number of time steps is reduced. Using large time slices is therefore not relevant as the final stage of the simulation falls into a low global energy level, which may be correctly handled by the coarse integrator itself. For large Reynolds numbers instead, a much larger number of time steps is required and Parareal can be applied with large time slices, while keeping a global energy level out of reach (regarding direct simulation constraints) for the coarse propagator.

					e k e	e
	80 46	1.4	4.4	26	16 -0.85% -2.8%
	160 124	2.7	4.1	41	30 -0.68% -2.0%
	320 322	3.4	4.1	55	46 -0.60% -1.5%
	ergy and on the dissipation are given in Tab. 3. The
	low values of both relative errors indicate that per-
	forming only one iteration of Parareal is also mean-
	ingful when the Reynolds number is increased. This
	relative independence to the Reynolds number is at-
	tributed to the nature of the coarse integrator. The
	spatial coarsening challenges the ability of

Table 4 :

 4 Influence of N ∆ t , G on the relative errors on k e and after one iteration of Parareal with spatial coarsening (Re λ 0 = 322, injection and 7 th -order interpolation as transfer operators).

	N ∆ t , Ĝ	e k	e
	10 -0.24% -3.3%
	20 -0.36% -2.7%
	40 -0.49% -3.4%
	Ĝ alone -0.63% -4.2%

Table 5 :

 5 Influence of C T ,1 and C I on the relative gain of the space-time parallelization σ S,T after one iteration of Parareal with spatial coarsening ( K = 1). Case (a) includes only C T ,1 , whereas case (b) considers both C T ,1 and C I in relation (32) (N = 4, N ∆ t , Ĝ = 20, injection and 7 th -order interpolation as transfer operators).

	N proc	Case σ S,T	σ S,T
	640	(a)	19%	24%
		(b)	18%	24%
	1280	(a)	47%	54%
		(b)	46%	54%

We define the parallel efficiency as the ratio of the parallel speedup over the number of processes involved.

Since there is no mean flow in this configuration, we cannot use the classical definition of the Reynolds number. We rather use Re λ (see e.g.[START_REF] Donzis | Dissipation and enstrophy in isotropic turbulence: Resolution effects and scaling in direct numerical simulations[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF]), which characterizes the ratio of turbulent eddy viscosity, over the fluid viscosity. This Reynolds number, at least one order of magnitude lower than the classical one, has therefore no relevance regarding laminar to turbulent transition, but instead informs about the range of length scales present in an already turbulent flow (see e.g.[START_REF] Pope | Turbulent Flows[END_REF] for a complete description

).[START_REF] Bal | Analysis of the convergence and of the stability of the Parareal algorithm for partial differential equations[END_REF] For isotropic turbulence λ is conveniently expressed as λ = u 15 µ /( ρ ) (see, e.g. ,[START_REF] Pope | Turbulent Flows[END_REF] Sec. 6.3, p.199]).

We refer the reader to only Sec. for numerical experiments using two iterations ( K = 2) of Parareal with spatial coarsening.

Acknowledgements The authors thank the two reviewers for their comments. The authors would like to acknowledge CALMIP for the dotation of computing hours on the EOS computer. This work was granted access to the HPC resources of CALMIP under allocation P1425 and P17005. Thibaut Lunet would like to thank Daniel Ruprecht, Michael Minion, Debasmita Samaddar, Jacob B. Schroder, Martin J.

Gander, Robert Speck, Andreas Schmidt, Felix Kwok, Beth Wingate and Raymond J. Spiteri for the interesting discussions at the 5th Workshop on Parallel-in-time integration at the Banff International Research Station (BIRS), and the Jülich Supercomputing Center for granting a Student Travel Award for attending this conference.

This research is partly funded by "Région Occitanie/ Pyrénées-Méditerranée" as part of the "Développement de nouvelles stratégies pour le calcul massivement parallèle à l'échelle exa en mécanique des fluides numérique (DENSE)" project.