J.-B Lagaert

G Balarac

G.-H Cottet

Hybrid spectral-particle method for the turbulent transport of a passive scalar

Keywords: Direct numerical simulation, Hybrid method, Particle method, Turbulent mixing, Passive scalar, Universal scaling

This paper describes a novel hybrid method, combining a spectral and a particle method, to simulate the turbulent transport of a passive scalar. The method is studied from the point of view of accuracy and numerical cost. It leads to a significative speed up over more conventional grid-based methods and allows to address challenging Schmidt numbers. In particular, theoretical predictions of universal scaling in forced homogeneous turbulence are recovered for a wide range of Schmidt numbers for large, intermediate and small scales of the scalar.

Introduction

The prediction of the dynamics of a scalar advected by a turbulent flow is an important challenge in many applications. The scalar field can be used to represent various quantities transported by the flow. In combustion, the mixture fraction is a conserved scalar used to describe mixing between fuel and oxidizer [START_REF] Pitsch | Large-eddy simulation of turbulent reacting flows[END_REF]. The prediction of scalar in environmental flows is also of great importance [START_REF] Michioka | High-Resolution Large-Eddy Simulations of Scalar Transport in Atmospheric Boundary Layer Flow over Complex Terrain[END_REF]. The temperature field is another type of advected scalar which is critical in many applications, e.g. to simulate the cooling systems used for nuclear reactors [START_REF] Kuczaj | Large-Eddy Simulation study of turbulent mixing in a T-junction[END_REF]. Passive scalars can finally be used to capture interfaces in multiphase flows [START_REF] Sethian | Level set methods for fluid interfaces[END_REF] or determine the dynamical properties of turbulent flows [START_REF] Conti | GPU and APU computations of finite time lyapunov exponent fields[END_REF].

A passive scalar, θ, is governed by an advection-diffusion equation,

∂θ ∂t + u • ∇θ = ∇ • κ ∇θ ( 1
)
where κ is the molecular scalar diffusivity and u the flow velocity field. The phenomenology of passive scalar convection-diffusion depends on the molecular Schmidt number, the viscosity-to-diffusivity ratio, Sc = ν/κ. For turbulent flow, the Kolmogorov scale, η K , is defined as the smallest length scale of the turbulent motion. Similarly, for Schmidt number higher than one, the Batchelor scale, η B , is defined as the smallest length scale of the scalar fluctuations. The Batchelor and Komolgorov scales are related by η B = η K / √ Sc. The Batchelor scale is thus smaller than the Kolmogorov scale. This means that, for Schmidt number larger than one, scalar dynamics can occur at scales smaller than the smallest velocity eddy, and therefore requires important computational resources. Donzis et al. [START_REF] Donzis | Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence[END_REF] performed DNS of turbulent transport by means of pseudo-spectral methods using up to 4096 modes in each direction to study universal scaling laws of a passive scalar.

In the above reference the same grid resolutions and numerical method are used for the momentum and the scalar equations. The two-scale nature of turbulent scalar transport, however, makes it natural to use different grids and different numerical methods for the scalar and the momentum. In a recent work, Gotoh et al. [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] describe a hybrid method combining a spectral method for the Navier-Stokes equation and compact finite-difference schemes for the scalar advection. This hybrid method is validated and applied for simulations of decaying turbulence at Schmidt numbers of 1 and 50. Significant gains were obtained in comparison with methods using spectral discretizations for both the momentum and the scalar.

For large Schmidt numbers, the scalar dynamics is essentially governed by advection, a regime for which Lagrangian or semi-Lagrangian methods are ideally suited. An important feature of these methods, which makes them particularly appealing in the case of high Schmidt numbers, is that they are stable under conditions that are related to the flow strain and not to the gridsize. In practice this means that the time-step used for the scalar equation will depend on the grid resolution used for the momentum equation even if a much finer grid is used for the scalar. Inspired by this observation, we proposed in [START_REF] Cottet | Subgrid particle resolution for the turbulent transport of a passive scalar[END_REF] to couple semi-Lagrangian particle methods at different gridresolutions for both the scalar transport and the Navier-Stokes equations. This reference provides a proof of concept that scalar spectra and structures are resolved with the same accuracy and much less computational effort in a hybrid method using a coarse resolution for the momentum than in a fully resolved high resolution method. This work was pursued in [START_REF] Lagaert | Particle method: an efficient tool for direct numerical simulations of a high Schmidt number passive scalar in turbulent flow[END_REF], to investigate the universal laws for large, intermediate and small scales of the scalar for Reynolds numbers (based on the Taylor micro-scale) between 80 and 160 and Schmidt numbers between 0.7 and 16. In this reference, a particle method for the scalar equation was coupled with a pseudo-spectral method for the Navier-Stokes equations.

The purpose of the present paper is to describe and validate the hybrid spectral-particle numerical approach used in [START_REF] Lagaert | Particle method: an efficient tool for direct numerical simulations of a high Schmidt number passive scalar in turbulent flow[END_REF], and to discuss its efficiency, in particular in comparison with fully resolved methods using spectral discretizations for both the scalar and momentum equations, and with the hybrid method proposed in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF].

An outline of this paper is as follows. In section 2, we describe the semi-Lagrangian particle method used for the scalar equation, the pseudo-spectral method used for the momentum equation and the coupling strategy. We also indicate the approach to run the hybrid method on massively parallel machines. In section 3, we test our method in decaying turbulence experiments similar to those in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] and discuss its accuracy, cost and overall efficiency. In section 4, we apply our method to investigate the physics of turbulent transport in forced homogeneous turbulence over a wide range of Schmidt numbers. Section 5 is devoted to concluding remarks and future directions that we are currently exploring.

Hybrid spectral-particle method

In this section we first describe the particle method used to solve the scalar equation, then the pseudo-spectral method used for the Navier-Stokes equation and the coupling strategy. We also explain our strategy to optimize the parallel performance of the hybrid method.

Semi-Lagrangian particle methods

The principle of particle methods for the advection of a given quantity is to concentrate this quantity on a set of particles and to follow these particles with the advection field. These methods are conservative by nature and free of CFL stability conditions. Continuous fields or grid values are recovered from the particles by mollification or interpolation [START_REF] Cottet | Vortex methods[END_REF]. The numerical analysis of these methods shows that a strong strain in the advection field can create distortions in the particle distribution and deteriorate the accuracy of the method. To overcome this difficulty, it is common practice to remesh particles on a regular grid through interpolation [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | Vortex methods[END_REF]. In the context of the advection of a vorticity field to solve the incompressible Navier-Stokes equation in vorticity form, these methods have been validated against spectral or finite-difference methods and applied in bluff body flows [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | Advances in direct numerical simulations of 3d wall-bounded flows by vortex-in-cell methods[END_REF][START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Application to the sphere at Re=300, 500 and 1000[END_REF][START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF], in homogeneous decaying turbulence [START_REF] Cottet | A comparison of spectral and vortex methods in three-dimensional incompressible flows[END_REF] and in vortex dynamics [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high reynolds numbers[END_REF][START_REF] Van Rees | Vortex tube reconnection at Re = 10 4[END_REF]. In the context of scalar advection they have been used for Lagrangian discretizations of level set methods [START_REF] Hieber | A lagrangian particle level set method[END_REF][START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF][START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] and for the determination of Lyapounov exponents of flow maps [START_REF] Conti | GPU and APU computations of finite time lyapunov exponent fields[END_REF].

When particles are remeshed at every time-step, which is often the case in practice, one obtains a class of conservative semi-Lagrangian methods that can be analyzed as CFL-free finite-difference methods [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF]. Remeshing is performed through interpolation. In one dimension it can be expressed by the following formula:

θ i = p θ p Λ( x i -x p ∆x θ ),
where Λ is the interpolation kernel, x i denote the grid points and x p the particle locations after advection. The summation concerns particles which belong to the support of the kernel around a given grid point. In the present paper, particles are advected by a second order Runge-Kutta scheme and we use the following kernel, derived in [START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF],

Λ(x) =          1 12 (1 -|x|) 25|x| 4 -38|x| 3 -3|x| 2 + 12|x| + 12 if 0 ≤ |x| < 1 24 (|x| -1) (|x| -2) 25|x| 3 -114|x| 2 + 153|x| -48 if 1 ≤ |x| < 1 24 (3 -|x|) 3 (5|x| -8) (|x| -2) if 2 ≤ |x| < 0 if 3 ≤ |x|.
(2) This kernel has a support of size 6. It is shown in [START_REF] Cottet | High-order semilagrangian particles for transport equation: numerical analysis and implementation issues[END_REF] that this remeshing kernel provides a scheme that is fourth order accurate in space for constant velocity field, and second order in the general case, provided the time-step is bounded by (max |∇u|) -1 . Throughout this paper we used the following value for the scalar time-step:

∆t θ = (max |∇ u|) -1 . (3)
Note that (max |∇ u|) -1 is the relevant time scale for the advection of a quantity with a velocity u. For scalar advection in three dimensions, following [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] we use a Strang splitting method, where particles are moved and remeshed alternatively in the 3 directions. More precisely, a complete advection-step with time-step ∆t θ of the scalar is performed by successively solving the advection equation along the x-axis, y-axis, z-axis, y-axis, x-axis, for ∆t θ /2, ∆t θ /2, ∆t θ , ∆t θ /2 and ∆t θ /2, respectively. After advection and remeshing of particles on the cartesian grid, scalar values are diffused on the grid. In the present study, we use an implicit spectral diffusion solver. In practice, with Schmidt numbers larger than one, we have observed that the diffusive time-scales were of the same order or bigger than the value given by (3). This time-step value is therefore consistent with the time-scales that need to be resolved to ensure accurate results for the advection-diffusion equations that we had to deal with.

Hybrid spectral-particle method

One of the advantages of semi-Lagrangian particle methods is that they are readily applicable to coupling approaches. In our case, we couple the advection-diffusion of a scalar by the semi-Lagrangian particle method just described with a pseudo-spectral method used to solve the Navier-Stokes equations. Both particle and spectral solvers are used in parallel, and employ different grid resolutions to take into account the separation between the Kolmogorov and Batchelor scales. The velocity field is computed by a pseudospectral solver with a second-order explicit Runge-Kutta time-advancement scheme, and a CFL number equal to 0.5. If ∆x u and ∆t u are the grid size and time-steps for Navier-Stokes solver, the CFL condition yields

∆t u = ∆x u 2 max |u| . (4)
The classical 3/2 rule is used for de-aliasing the non-linear inertial term of the Navier-Stokes equation and the viscous terms are treated exactly by integrating analytically the spectral form of the diffusion equation. The velocity obtained from the Navier-Stokes solver is interpolated in spectral space to obtain particle velocities. Although less expensive interpolations could have been used, spectral interpolation has been chosen because it minimizes interpolation errors. Since the velocity resolution is taken in a DNS range, the spectral interpolation beyond the Kolmogorov scale is indeed nearly exact. Moreover our numerical experiments showed that the computational overhead resulting from this interpolation method remained small in comparison with the computational cost associated with the scalar and momentum equations.

In practice, due to (3) and ( 4), the time-step used for the scalar equation is always larger than the time-step used for the Navier-Stokes equation. ∆t θ can be chosen to be an integer multiple of ∆t u and several sub-steps of the Navier-Stokes solver are performed for one iteration of the scalar equation.

Parallel efficient implementation

Our goal is to perform direct numerical simulations (DNS) of turbulent transport at high Schmidt numbers. A very fine spatial resolution is therefore required, which makes the development of algorithms optimized for massively parallel computers mandatory.

Parallel algorithm for the pseudo-spectral solver

The pseudo-spectral solver uses a classical parallel design based on twodimensional subdivisons. The 3D fast Fourier transform is split in three one-dimensional fast Fourier transforms. The parallelization is done by subdividing the 3D geometry along only 2 directions: along Y and Z in the real space and along X and Y on the spectral space (see [START_REF] Lagaert | Particle method: an efficient tool for direct numerical simulations of a high Schmidt number passive scalar in turbulent flow[END_REF] for more details). In order to achieve the coupling between the two solvers, the domain subdivision used for the particle method implementation is the same as the one used inside the pseudo-spectral library in the real space. This avoids data-reorganization.

Parallel design of the particle method

Our strategy to limit communications in the particle method is to define work items as a group of particles, defined on the basis of their location and a geometric subdivision, that are followed during a time-step (see algorithm1). The directional splitting further dictates the choice of a (large) group of particle lines for these work items, which allows to gather communications (lines 8, 9 and 14 of algorithm 1). Each time step requires only four communications (one at lines 8 and 9 and 2 at line 14 of algorithm 1) per line group for each MPI-process.

Overall scalability

The parallel scalability has been tested on different architectures. Table 1 presents results on a Blue-Gene/Q cluster and on an Intel-based cluster using quad-cores with Intel Westmere-EP processors. Due to memory limitations, the grid sizes are not the same on the two clusters. For the purpose of these scalability measures, the Navier-Stokes and scalar equations are solved with Algorithm 1 Advection of a particle line.

1: for each (1D-) line do 2:

Create a particle on each mesh point and initialize its position:

p i n = i • ∆x θ . 3:
Compute the intermediate RK2 positions for each particle:

p i n+1/2 = p i n + ∆t θ /2 • v i n .

4:

Compute the integer truncation I i n of p i n+1/2 /∆x.

5:

for m = 0 to 1 do

6:

Compute the rank R m (i, n) of the process which stored the velocity on node

I i n + m. 7: if my MPI-rank is different from R m (i, n) then 8:
Send a MPI message to process R m to ask for the velocity on node I i n + m.

9:

Receive the velocity on node I i n + m from process R m (i, n). Interpolate the velocity on position p i n+1/2 from velocity on node I i n and (

I i n + 1) {as p i n+1/2 ∈ [I i n ; I i n + 1]}.
13:

Compute the final particle positions:

p i n+1 = p i n + ∆t θ /2 • v n (p i n+1/2 ).
14:

Remesh particles. Table 1: Strong scalability of the hybrid method on different architectures.

To compute the scalability, the run on the minimal cores number is used as a reference. N u and N θ denote the number of mesh points for velocity and scalar, respectively.

Validation and performances of the hybrid method

Simulation setup

To validate and to compare the hybrid spectral-particle method with, on the one hand, a pure spectral method and, on the other hand, the hybrid method proposed by Gotoh et al. [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF], we have implemented the same testcases as in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF]. In all cases, we solved the Navier-Stokes equations with a classical pseudo-spectral method using a second order Runge-Kutta timestepping and the CFL condition [START_REF] Sethian | Level set methods for fluid interfaces[END_REF].

In these test-cases, both scalar and velocity freely decay from random initial fields. The random fields are multivariate Gaussians with a given spectrum. Following [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] these initial spectra, for the velocity and scalar, are chosen as follows

E(k, t = 0) = 16 2 π k -5 0 k 4 exp (-2(k/k ν ) 2 ), (5)
E θ (k, t = 0) = 32 2 π k -5 0 k 4 exp (-2(k/k θ ) 2 ), (6)
with k ν = k θ = 6 and k 0 = 5. In this section we report simulations corresponding to two different Schmidt numbers, 1 and 50. As in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF], the value of the viscosity ν was adjusted to keep the same scalar resolutions in both cases.

Table 2 presents the physical parameters and the grid resolutions used in our different simulations, and table 3 provides the statistics that determine the initial conditions. Here the Reynolds and Péclet number are defined from the Taylor micro scales,

R λ = u λ ν λ = u 2 1 (∂u 1 /∂x 1 ) 2 , (7)
P λ = u λ θ κ λ θ = θ 2 (∂θ/∂x 1 ) 2 . ( 8
)
For completeness, we recall the definition of the mean kinetic energy and scalar variance dissipation rates,

¯ = ν 2 i,j ∂u j ∂x i + ∂u i ∂x j 2 , (9) χ
= 2κ i ∂θ ∂x i 2 . ( 10
)
We will also consider the mixed skewness S uθ of the velocity and scalar gradients, which measures the transfer of the scalar excitation toward small scale

S uθ = ∂u 1 /∂x 1 (∂θ/∂x 1 ) 2 (∂u 1 /∂x 1 ) 2 1/2 (∂θ/∂x 1 ) 2 , ( 11
)
and the flatness of the scalar gradient

F ∂θ/∂x = (∂θ/∂x 1 ) 4 (∂θ/∂x 1 ) 2 2 . ( 12
)

Accuracy of the hybrid method

In all the following comparisons, the time is normalized by the eddy turnover time. As it is customary, the spatial resolution of our DNS is chosen such that K u max η K and K θ max η B larger than 1.5 [START_REF] Pope | Turbulent Flows[END_REF], where K u max and K θ max are the maximal wave number for the discrete velocity and the scalar, respectively.

Run

k E(k, t) Reference -N u = N θ = 1536 3 Spectral -N u = N θ = 1024 3

Test case I: Sc = 1

For Sc = 1, since η K ≈ η B , the same resolution is needed for velocity and scalar fields. A reference spectral simulation has been run with 1536 3 points. Figure 1 shows that a spectral simulation with 1024 3 points already provides a well resolved DNS, at least until time t = 4. This is confirmed by figure 2, which shows the kinetic energy spectrum at several times for both resolutions.

To study the accuracy of the hybrid scheme, the spectral simulations are now compared with simulation using the hybrid spectral-particle method on the same mesh (N θ = 1024 3 points) and with finer resolutions using N θ = 1280 3 = (1.25 × 1024) 3 or N θ = 1536 3 = (1.5 × 1024) 3 points for the scalar (see Table 2).

Figure 3 shows first the scalar variance spectrum for the different simulations, at various time. All the spectra well agree except for run 3 (hybrid method with N θ = 1024 3 ) at the smallest scales. This means that the particle discretization of the scalar equation leads to a slight numerical dissipation at the smallest scales. Note that this has also been observed in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] with a fourth order finite difference scheme. Figure 3 also shows that when K θ max η B = 2.25, which corresponds to 1.5 times the classical resolution criterion, the hybrid method recovers most of the spectrum.

Additional comparisons can be performed based on the probability density functions (PDF). gradient scalar and of the scalar dissipation rate, respectively. The probability distribution functions of the scalar produced in runs 2 and 3 are in excellent agreement. The PDFs of the scalar gradient and of the scalar dissipation are more discriminant, and some discrepancies with the reference solution are observed near the tails of the distributions even for the spectral simulation using 1024 3 grid points. For the hybrid simulations, the highest resolution (N θ = 1536 3 ) well reproduces the reference results, whereas the smaller resolutions (N θ = 1024 3 and 1280 3 ) lead to a slight under-prediction of the tails, comparable to what is observed on the 1024 3 spectral simulation. The PDF of the mean dissipation follows a exp(-c(χ/ χ) α ) law [START_REF] Chertkov | Intermittent dissipation of a passive scalar in turbulence[END_REF].

The exponent α is expected to slowly decay in time, but it is interesting to evaluate its value and compare it with previous studies. From our results, we obtain α = 0.354, 0.370, 0.372, 0.365 and 0.357 at u 0 k 0 t = 6 for run 1 to 5, respectively. These values are close to the theoretical prediction α = 1/3 of the Kraichnan model [START_REF] Chertkov | Intermittent dissipation of a passive scalar in turbulence[END_REF]. Some complementary statistics, similar to those presented in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF], are finally presented. The time evolution of the scalar mean dissipation [START_REF] Cottet | Vortex methods[END_REF] and of the flatness (12) of the scalar gradient (fig. 7 and 8 respectively) exhibit some discrepancies between the spectral method and the hybrid method using the same scalar resolution. For N θ ≥ 1280 3 (run 4 and 5), the hybrid method provides a very good agreement with the spectral method. The last statistic is the mixed skewness of velocity and scalar gradient S uθ [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. Its time evolution is shown by figure 9. As explained in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF], the reference run slightly over predicts this diagnostics: in all cases, the finite resolution limits the energy transfer near the cut-off. This resolution effect makes it harder to compare the different results, but it appears that the hybrid method delivers very accurate results for N θ ≥ 1280 3 . From all these data, on can conclude that, for S c = 1, the hybrid spectralparticle method with N θ = 1536 3 (run 5), although slightly over-dissipative near the cut-off wave number, exhibits for all diagnostics an accuracy similar to the pseudo-spectral method and to the high order compact finitedifference methods proposed in [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] . In other words, the modified DNS criterion, K θ max η B ≥ 2.25 = 1.5 2 ensures that the hybrid method can be run in a DNS regime. For 1.5 < K θ max η B < 2.25, run 4 shows that the statistic quantities are rather well recovered. We will see below (section 3.3) that, even for this small value of the Schmidt number and with higher resolution, the hybrid method leads to a computational speed-up over both the purely spectral method and the hybrid method of [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF].

Test-case II: Sc = 50

To address this case, following [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] the fluid viscosity is decreased to 1.2 × 10 -2 . This allows to respect the DNS criterion with a resolution of 1024 3 points for the scalar and 256 3 points for the velocity (see Fig. 10). For the hybrid method, a simulation with 1280 3 points for the scalar (run 7) is performed, based on the modified DNS criterion K θ max η B = 2.25. Another simulation (run 6) is also performed with 1024 3 points to compare spectral and hybrid methods on the same resolution (see table 2).

u 0 k 0 t K u max η K , K θ max η B Spectral -K u max η K -N u = 256 3 Spectral -K θ max η B -N θ = 1024 3
First, figure 11 shows the scalar variance spectra at time u 0 k 0 t = 1.6 and 3. The hybrid method using N θ = 1280 3 particles is in excellent agreement with the spectral method except near the cut-off, whereas, as expected, the run 6 shows some numerical discrepancy at high wave numbers. The PDFs of the scalar gradient and the scalar dissipation are shown on Figures 12 and13. All the runs provide the same scalar gradient distribution. The hybrid method accurately reproduces the tails of the PDF of the scalar dissipation rate only for N θ ≥ 1280 3 . The fit with the theoretical form exp c(χ/ χ) α gives α = 0.351, 0.374 and 0.350 at time k 0 u 0 t = 3 for the spectral method and the hybrid method with N θ = 1024 3 and 1280 3 , respectively. Finally, figures 14 and 15 show the mixed skewness and the flatness of the scalar gradient. The results obtained with the spectral method and the hybrid method with N θ = 1280 (i.e. K θ max η B ≥ 2.25) are in excellent agreement. These experiments confirm that the condition K θ max η B ≥ 2.25 on the scalar resolution (coupled with the condition K u max η K ≥ 1.5 to ensure that the Navier-Stokes equations are solved in a DNS regime) ensures accurate DNS with the hybrid method. For 1.5 < K θ max η B < 2.25, the global statistic quantities are correct but the smallest scales are slightly underestimated. Nevertheless we will see in the next section that this error remains small and does not prevent to study qualitatively the spectrum decay and discriminate between universal scaling laws. We now explain how this hybrid strategy pro- (∂θ/∂x 1 )/σ σP(∂θ/∂x 1 )

Figure 12: PDF of the scalar gradient at time u 0 k 0 t = 3 for Sc = 50. σ is the mean square root of the gradient scalar (computed in spectral run). See figure 11 for legend. vides important computational savings even when a finer spatial resolution is used for the scalar.

Numerical efficiency

Table 4 presents the computational cost of the different methods. In this table, the total CPU time to reach a given simulation time u 0 k 0 t = 6 is reported. We recall that, in practice, several sub-steps of the Navier-Stokes solver are performed for one iteration of the scalar equation, since the particle method can be used with large time-steps. To properly evaluate the computational cost of the particle method and the computational efficiency resulting from the ability to use large time-steps, additional runs (3bis, 4bis, ...) have been performed with the hybrid method using the same (small) time-steps as the Navier-Stokes equations. For these runs, the CPU time per time-step is also reported.

All the simulations reported in Table 4 have been performed on 2048 cores of a IBM Blue Gene Q. Note that this computer exhibits a rather low CPU frequency [START_REF] Pitsch | Large-eddy simulation of turbulent reacting flows[END_REF][START_REF] Donzis | Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence[END_REF] GHz for each core). Table 4: Numerical efficiency of the different methods -Runs are performed on 2048 cores of a Blue Gene Q. N u , N θ denotes the spatial resolution for velocity and scalar and ∆t u , ∆t θ are the numerical time steps for momentum and scalar equations.

A first conclusion is that the particle method for the scalar equation runs about 10% faster than the spectral method using the same grid resolution and time-sep (runs 2 and 3bis). On Intel-based clusters, which are more common, similar comparisons have been performed, and the acceleration factor reaches the value of 2. However the most important speed-up factor comes from the large time step allowed by the hybrid method. As expected, the scalar time step does not depend on the scalar spatial resolution. For test-cases at Sc = 1, the same time-step is used for N θ = 1024 3 , 1280 3 and 1536 3 whereas a spectral method, or any grid-based method, would require to divide it by 1.25 for N θ = 1280 3 and 1.5 for N θ = 1536 3 . Note that in runs 3, 4 and 5, the scalar time step is not constant as the maximum value of the velocity gradient varies during the simulation: at the beginning, the velocity gradient increases until the kinetic energy spectrum reaches the smallest scales (near the Kolmogorov scale) and then it starts to decrease due to the energy dissipation. The range of the variations of ∆t θ is presented in table 4. ∆t θ always remains much larger than ∆t u . Comparisons of the total CPU times for runs 2, 3, 4 and 5 show that, even when the spectral method is used under the condition K θ max η B ≈ 1.5, and the spectral-particle scheme under the more severe condition K θ max η B ≥ 2.25, the hybrid method is about 20% faster.

The efficiency of the hybrid method is even more pronounced for larger values of the Schmidt number. For Sc = 50, in hybrid methods based on Eulerian schemes [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF] the time-step for the scalar equations will be determined by the scalar resolution. This is not the case for the spectral-particle hybrid method. A comparison of the total CPU time of runs 6 and 9, shows that, even with a higher resolution, the hybrid spectral-particle method leads to a speed up factor of 32 over the purely spectral method. These results can be compared to those obtained by the hybrid spectral-finite-difference method of [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF]. This reference reports a speed up factor about 4 for hybrid spectralfinite-difference methods in comparison with a purely spectral method. Although this comparison must be taken with care, due to the difference in the machines and of the spectral codes, it demonstrates the value of the hybrid spectral particle method for high values of the Schmidt number.

Application to the advection of passive scalar at high Schmidt numbers

The hybrid spectral-particle method is now used to study the turbulent mixing on a wide range of Schmidt numbers, with numerical parameters that have been validated in the previous section. Theoretical studies explain the influence of the Schmidt number on the behavior of the scalar variance spectrum [START_REF] Lesieur | Turbulence in fluids, Fluid mechanics and its applications[END_REF]. For a Schmidt number larger than one, Batchelor [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF] reports that the classical Corrsin-Obukhov cascade associated with a k -5/3 law (where k is the wave number) for the scalar variance spectrum [START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF][START_REF] Obukhov | The structure of the temperature field in a turbulent flow[END_REF] is followed by a viscous-convective range with a k -1 power law. This viscous-convective range is followed by the dissipation range, where various theoretical scalings have been proposed for the spectrum [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF][START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF]. The goal of this section is to compare DNS results based on the hybrid-spectral method with these theoretical predictions for small, intermediate and large scales and for various Schmidt numbers.

Flow configuration

Several simulations have been performed in the context of forced homogenous isotropic turbulence, in a 3D periodic box with a length 2π. The forcing scheme used to obtain a statistical steady flow follows the one proposed by Alvelius [START_REF] Alvelius | Random forcing of three-dimensional homogeneous turbulence[END_REF]. To achieve a steady state for the scalar, a forcing scheme is also applied to low wave number modes in Fourier space, similarly to the velocity Table 5: Setup of simulations performed. ∆t u is the time step used to solve the Navier-Stokes equation with a pseudo-spectral solver. ∆t θ is the time step used to solve the scalar transport equation with the particle method. ∆t θ spec is the time step which would be needed if a pseudo-spectral method was used for the same number of scalar grid points.

R λ N u K u max η K ∆t u Sc N θ K θ max η B ∆t θ ∆t θ
forcing [START_REF] Silva | Analysis of the gradient-diffusion hypothesis in large-eddy simulations based on transport equations[END_REF]. Two Taylor-scale Reynolds numbers, R λ , are considered, 130 and 210, using a resolution of N u = 256 3 and 512 3 grid points, respectively. Figure 16 (left) shows the compensated spectrum for the kinetic energy, where ¯ is the mean energy dissipation rate. For the scalar field, the mesh resolution is increased with the Schmidt number [START_REF] Donzis | Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence[END_REF] from 512 3 to 3064 3 grid points. Simulation details are given in Table 5. The time-steps used in these simulations are independent of the scalar resolutions, which allowed important cost reduction for the high Schmidt number cases. On the basis of the comparisons reported in the previous section, given the values of the time-steps shown in Table 5, for Sc = 128 the speed-up provided by the spectral-particle method over a purely spectral method can be estimated at about 100. Figure 17 illustrates the scales separation between the Kolmogorov and Batchelor scales for the highest Schmidt number case, Sc = 128.

Scalar spectrum analysis

The behaviors of the scalar variance spectrum are next studied at large, intermediate and small scales from this DNS database and compared with theoretical predictions. First, at the scales beyond the forcing peak, the classical Corrsin-Obukhov cascade is expected to characterize the inertialconvective range. Similarly to the inertial range of the kinetic energy spec- trum, it is expected that this cascade follows a k -5/3 law [START_REF] Obukhov | The structure of the temperature field in a turbulent flow[END_REF][START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF]. Figure 16 (right) shows the scalar spectra for various Schmidt numbers, compensated by the Corrsin-Obukhov scaling, where χ is the mean scalar dissipation rate. As expected, the results show an inertial-convective range independent of the Schmidt number. The k -5/3 range of the scalar spectrum is found more clearly than the k -5/3 range of the energy spectrum. Indeed, the scalar spectrum exponent is known to tend to the -5/3 value more rapidly than the energy spectrum exponent [START_REF] Lee | Scaling range of velocity and passive scalar spectra in grid turbulence[END_REF][START_REF] Danaila | Spectrum of a passive scalar in moderate Reynolds number homogeneous isotropic turbulence[END_REF]. Note that the end of the inertialconvective range appears around kη K ≈ 0.1 (which is roughly the Taylor scale) independently of the Schmidt number.

Beyond this range, for Schmidt numbers larger than one, Batchelor [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF] described the development of the viscous-convective range with a k -1 law. This scaling is due to the velocity small scales strain effect on the scalar field. Figure 18 shows the scalar spectra for various Schmidt numbers, compensated by the Batchelor's scaling. For clarity, the results are shown with wavenumbers multiplied by the Kolmogorov scale (Figure 18, left picture) or by the Batchelor scale (Figure 18, right picture). From our numerical results, we observe that the k -1 power law starts from kη K ≈ 0.1, after the inertial-convective range, and this viscous-convective range grows with the Schmidt number (Figure 18, left picture).

The form of the scalar variance spectrum in the dissipation range (following the viscous-convective range) is also studied (Figure 18, right picture). Two distinct theoretical behaviors have been proposed by Batchelor [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF] and Kraichnan [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF]. For scales beyond the Kolmogorov scale, Batchelor [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF] assumes that the scalar at high Schmidt number is strained by the smallest scales of the velocity. Considering spatial fluctuations of the local strain rates, but neglecting temporal fluctuations, Batchelor predicts the following form of the scalar variance spectrum

E θ (kη B ) χ ν/ η B = q kη B exp -q(kη B ) 2 , ( 13
)
where q is a constant. Alternatively, Kraichnan [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF] assumes that strain rates rapidly fluctuate in time and obtains

E θ (kη B ) χ ν/ η B = q kη B 1 + (6q) 1/2 kη B exp -(6q) 1/2 kη B . (14)
Note that the distinction between these behaviors has some practical implications on our understanding of energy transfer between ocean and atmosphere [START_REF] Bogucki | Experimental evidence of the kraichnan scalar spectrum at high reynolds numbers[END_REF], for example. Although the different curves do not perfectly coincide, our DNS results clearly show a good agreement with the Kraichnan form, even for the runs where 1.75 < K θ max η B < 2.25. This result confirms previous studies performed in [START_REF] Bogucki | Experimental evidence of the kraichnan scalar spectrum at high reynolds numbers[END_REF][START_REF] Donzis | The batchelor spectrum for mixing of passive scalars in isotropic turbulence[END_REF][START_REF] Bogucki | Direct numerical simulations of passive scalars with P r > 1 advected by turbulent flow[END_REF]. Note that, Donzis et al. [START_REF] Donzis | The batchelor spectrum for mixing of passive scalars in isotropic turbulence[END_REF] observed that, even for Schmidt number smaller than 1, the scalar spectra collapse in the dissipative range when they are normalized by Batchelor variables. The same trend is found with our numerical results for Sc = 0.7 (Fig. 18 right). However, for Sc = 0.7 (the smallest Schmidt number value shown on figures), the spectrum does not exhibit a k -1 range (Fig. 18, left). This range is only found for Sc > 1, as predicted by theory.

Effect of forcing schemes

Our numerical results find the same behavior for the smallest scales as previous studies and observations [START_REF] Bogucki | Direct numerical simulations of passive scalars with P r > 1 advected by turbulent flow[END_REF][START_REF] Yeung | Schmidt number effects on turbulent transport with uniform mean scalar gradient[END_REF]. However, our numerical results also clearly show a -5/3 inertial-convective range for moderate Reynolds numbers, in contrast with the results reported in [START_REF] Donzis | The batchelor spectrum for mixing of passive scalars in isotropic turbulence[END_REF], where an imposed mean scalar gradient was used [START_REF] Yeung | Schmidt number effects on turbulent transport with uniform mean scalar gradient[END_REF]. To better understand the influence of the forcing scheme at large and intermediate scales, additional simulations have been performed with an imposed mean scalar gradient to maintain the scalar variance. Figure 19 compares the scalar variance spectra for two Reynolds numbers, 130 and 210 and two Schmidt numbers, 0.7 and 4. The spectra are compensated by using the scaling of the inertial-convective range. As expected the small scales behavior is not influenced by the forcing schemes. But, for moderate Reynolds number, the influence of the forcing schemes clearly appears. The simulations using a mean scalar gradient forcing have no clear k -5/3 range, and they have a large k -1 range beginning at large scales. In particular, for Sc = 0.7, the scalar spectrum exhibits a viscousconvective range for the mean scalar gradient forcing, in contrast with the results obtained with a low wave numbers forcing and with the theoretical prediction. When the Reynolds number increases, an inertial-convective range begins to appear for simulations with mean scalar gradient forcing, as shown by the spectra for R λ = 210. Note that Watanabe and Gotoh [START_REF] Watanabe | Intermittency in passive scalar turbulence under uniform mean scalar gradient[END_REF] have already observed an influence of the forcing scheme on the scaling exponents of the structure function.

Conclusion

In this paper we have presented a hybrid scheme for the turbulent transport of a passive scalar, combining a spectral method for the Navier-Stokes equations and a semi-Lagrangian particle method for the scalar transport. The method takes advantage of the Lagrangian nature of particle methods to enable high resolution of the scalar at affordable numerical cost.

A validation study in decaying homogeneous turbulence against a purely spectral method allowed to clarify the particle resolution necessary to resolve the smallest scales. This study showed the gain of the present method over hybrid methods combining spectral and finite-difference methods in the case of high Schmidt numbers. The efficiency of the method was further demonstrated in a study of universal scalar laws in forced homogeneous turbulence for Schmidt numbers up to 128. Using scalar resolution up to 3064 3 particles on massively parallel clusters, the method allowed to recover the theoretical predictions at large, intermediate and small scales.

Further works are currently undertaken to extend the capabilities of the method in several directions. The local nature of particle methods allows us to envision hybrid method where the scalar is solved on GPU along the line of [START_REF] Cottet | High-order semilagrangian particles for transport equation: numerical analysis and implementation issues[END_REF]. Moreover, thanks to the dimensional splitting, using high order particle remeshing schemes derived in [START_REF] Cottet | High-order semilagrangian particles for transport equation: numerical analysis and implementation issues[END_REF] should further improve the accuracy of the hybrid method for a marginal additional cost. Finally, the present study opens the way to a multi-scale approach to vorticity transport as an LES tool for the Navier-Stokes equations.

Figure 1 :

 1 Figure 1: Time evolution of K u max η K for Sc = 1.

Figure 2 :

 2 Figure 2: Kinetic energy spectra E(k, t) at time u 0 k 0 t = 0, 1.5, 3 and 6.

  Figures 4,[START_REF] Conti | GPU and APU computations of finite time lyapunov exponent fields[END_REF] and 6 show PDFs of the scalar, of the -N u = 1024 3 , N θ = 1024 3 Hybrid -N u = 1024 3 , N θ = 1280 3 Hybrid -N u = 1024 3 , N θ = 1536

Figure 3 :

 3 Figure 3: Spectra of the scalar variance E θ (k, t) at times u 0 k 0 t = 1.5 and 6.

Figure 4 :

 4 Figure4: PDF of the scalar value at time u 0 k 0 t = 3 and u 0 k 0 t = 6 for Sc = 1. σ denotes the mean square root of the scalar (computed in run 1). See figure3for legend.

Figure 5 :Figure 6 :

 56 Figure5: PDF of the scalar gradient at time u 0 k 0 t = 6 for Sc = 1. σ is the mean square root of the scalar gradient (computed in run 1). See figure3for legend.

Figure 7 :

 7 Figure 7: Time evolution of the mean scalar dissipation rate for Sc = 1. See figure 3 for legend.

Figure 8 :Figure 9 :

 89 Figure 8: Time evolution of the flatness of the scalar gradient for Sc = 1. See figure 3 for legend.

Figure 10 :

 10 Figure 10: Time evolution of K u max η K and K θ max η B for Sc = 50.

Figure 11 :

 11 Figure 11: Spectra of the scalar variance E θ (k, t) at time u 0 k 0 t = 1.6 and 3.0 for Sc = 50.

Figure 13 :Figure 14 :

 1314 Figure 13: PDF of the scalar dissipation at time u 0 k 0 t = 3 for Sc = 50. See Fig. 11 for legend.

Figure 15 :

 15 Figure 15: Time evolution of the flatness of the scalar gradient for Sc = 50. See Fig. 11 for legend.

Figure 16 :

 16 Figure 16: Left picture: Energy spectra for both Reynolds numbers, R λ ≈ 130 and R λ ≈ 210. The spectra are compensated by the Kolmogorov scaling. The arrow shows the direction of increasing Reynolds numbers. Right picture: Scalar variance spectra for R λ ≈ 130. The spectra are compensated by the scaling proposed by Corrsin-Obukov. The arrow shows the direction of increasing Schmidt numbers.

Figure 17 :

 17 Figure 17: x-y plan colored by the vorticity magnitude (left, blue regions are for the lowest vorticity values and red regions are for the highest vorticity values) and by the passive scalar (middle, blue regions are for the lowest scalar values and red regions are for the highest scalar values) for R λ ≈ 130 and Sc = 128. The zooms (right) correspond to the white box with a length of 8η K for the vorticity magnitude (top) and the scalar (bottom).

Figure 18 :

 18 Figure 18: Scalar variance spectra for R λ ≈ 130. The spectra are compensated by the scaling proposed by Batchelor for Schmidt number higher than one. The arrow shows the direction of increasing Schmidt numbers. For a sake of clarity, the results are shown with wavenumbers multiplied by the Kolmogorov scale (left) and by the Batchelor scale (right). For the dissipative region, the circles show the law proposed by Kraichnan and the squares show the law proposed by Batchelor (right).

Figure 19 :

 19 Figure 19: Scalar variance spectra for two Schmidt number, Sc = 0.7 and Sc = 4 for R λ ≈ 130 (left) and R λ ≈ 210 (right). The spectra are compensated by the scaling proposed by Corrsin-Obukov. The arrow shows the direction of increasing Schmidt numbers. The dashed and solid lines show the spectra for the low wave numbers forcing and constant gradient forcing, respectively. The dashed-dotted line corresponds to the k -1 scaling under the chosen normalization.

  15: end for the same time-step. Grid resolutions for the velocity and scalar used in this study are typical of those used for high Schmidt number simulations. Strong scalability on IBM Blue Gene Q for N θ = 2048 3 and N u = 256 3 . Strong scalability on Intel cluster for N θ = 1024 3 and N u = 128 3 .

	Cores 1024 2048 4096 8192 (a) Cores Time/iteration Scalability 102.02s -52.39s 97% 26.80s 95% 14.32s 89% 128 256 512 (b)	Time/iteration 55s 27.5s 14.4s	Scalability -100% 95%

Table 2 :

 2 Numerical parameters for the different runs.

		Scheme		ν	Sc	N u	N θ
	1 2 3 4 5	Reference (spectral) Spectral Hybrid Hybrid Hybrid	4.0 × 10 -4 4.0 × 10 -4 4.0 × 10 -4 4.0 × 10 -4 4.0 × 10 -4	1 1 1 1 1	1536 3 1024 3 1024 3 1024 3 1024 3
	6 7 8 9	Spectral Spectral Hybrid Hybrid		1.2 × 10 -2 1.2 × 10 -2 1.2 × 10 -2 1.2 × 10 -2	50 50 50 50	1024 3 256 3 256 3 256 3
		Sc	¯	χ	R λ	P λ
	Run 1-5	1	0.06498	0.2596	881	690
	Run 6-9	50	1.950	0.1557	29	1150

Table 3 :

 3 Statistics at initial time.

Acknowledgments

This work was partially supported by the Agence Nationale pour la Recherche (ANR) under Contracts No. ANR-2010-JCJC-091601 and ANR-2010-COSI-0009. G.-H. C. is also grateful for the support from Institut Universitaire de France. This work was performed using HPC resources from GENCI-IDRIS (Grant 2012-020611). The authors thankfully acknowledge the hospitality of the Center for Turbulence Research, NASA-Ames and Stanford University, where part of this work has been done during the Summer Program 2012. They are also grateful to J. A. Domaradzki for fruitful discussions.