Cristian C Lalescu
email: cristian.lalescu@mpcdf.mpg.de

Bérenger Bramas

Markus Rampp

Michael Wilczek

An Efficient Particle Tracking Algorithm for Large-Scale Parallel Pseudo-Spectral Simulations of Turbulence

Keywords: Particle tracking, Turbulence simulations, MPI, OpenMP

Particle tracking in large-scale numerical simulations of turbulent flows presents one of the major bottlenecks in parallel performance and scaling efficiency. Here, we describe a particle tracking algorithm for large-scale parallel pseudo-spectral simulations of turbulence which scales well up to billions of tracer particles on modern high-performance computing architectures. We summarize the standard parallel methods used to solve the fluid equations in our hybrid MPI/OpenMP implementation. As the main focus, we describe the implementation of the particle tracking algorithm and document its computational performance. To address the extensive inter-process communication required by particle tracking, we introduce a task-based approach to overlap point-to-point communications with computations, thereby enabling improved resource utilization. We characterize the computational cost as a function of the number of particles tracked and compare it with the flow field computation, showing that the cost of particle tracking is very small for typical applications.

Introduction

Understanding particle transport in turbulent flows is fundamental to the problem of turbulent mixing [START_REF] Taylor | Diffusion by continuous movements[END_REF][START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF][START_REF] Yeung | Lagrangian statistics from direct numerical simulations of isotropic turbulence[END_REF][START_REF] Yeung | Lagrangian investigations of turbulence[END_REF][START_REF] Toschi | Lagrangian Properties of Particles in Turbulence[END_REF][START_REF] Homann | Clustering of passive impurities in magnetohydrodynamic turbulence[END_REF] and relevant for a wide range of applications such as dispersion of particles in the environment [START_REF] Stohl | A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements[END_REF][START_REF] Stohl | Characteristics of atmospheric transport into the arctic troposphere[END_REF][START_REF] Behrens | Model simulations on the long-term dispersal of 137 Cs released into the Pacific Ocean off Fukushima[END_REF][START_REF] Haszpra | Dispersion of aerosol particles in the free atmosphere using ensemble forecasts[END_REF], the growth of cloud droplets through collisions [START_REF] Shaw | Particle-turbulence interactions in atmospheric clouds[END_REF][START_REF] Bodenschatz | Can we understand clouds without turbulence?[END_REF][START_REF] Devenish | Droplet growth in warm turbulent clouds[END_REF][START_REF] Grabowski | Growth of cloud droplets in a turbulent environment[END_REF][START_REF] Pumir | Collisional aggregation due to turbulence[END_REF], and phytoplankton swimming in the ocean [START_REF] Durham | Turbulence drives microscale patches of motile phytoplankton[END_REF][START_REF] Breier | Emergence of phytoplankton patchiness at small scales in mild turbulence[END_REF][START_REF] Pujara | Rotations and accumulation of ellipsoidal microswimmers in isotropic turbulence[END_REF]. Direct numerical simulations (DNS) of turbulence are nowadays an established tool for investigating such phenomena and have a long history in scientific computing [START_REF] Orszag | Numerical simulation of three-dimensional homogeneous isotropic turbulence[END_REF][START_REF] Yeung | An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence[END_REF][START_REF] Yokokawa | 4-tflops direct numerical simulation of turbulence by a Fourier spectral method on the earth simulator[END_REF][START_REF] Ishihara | Energy spectrum in high-resolution direct numerical simulations of turbulence[END_REF][START_REF] Buaria | A highly scalable particle tracking algorithm using partitioned global address space (pgas) programming for extreme-scale turbulence simulations[END_REF]. DNS have become a major application and technology driver in highperformance computing since the scale separation between the largest and the smallest scales increases drastically with the Reynolds number R λ , which characterizes the degree of smallscale turbulence [START_REF] Pope | Turbulent Flows[END_REF]. Dimensional estimates of the required computational resources scale at least as R 6 λ [START_REF] Pope | Turbulent Flows[END_REF]. Recent literature [START_REF] Yeung | Effects of finite spatial and temporal resolution in direct numerical simulations of incompressible isotropic turbulence[END_REF], however, shows that, due the occurrence of extremely small-scale structures, resolution requirements increase even faster than simple dimensional arguments suggest. Until today DNS have reached values of up to R λ ≈ 2300 [START_REF] Ishihara | Smallscale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient statistics[END_REF][START_REF] Yeung | Extreme events in computational turbulence[END_REF][START_REF] Ishihara | Energy spectrum in high-resolution direct numerical simulations of turbulence[END_REF], still smaller than the latest experiments, which have reached R λ > 5000 [START_REF] Küchler | Experimental study of the bottleneck in fully developed turbulence[END_REF], or natural settings such as cumulus clouds, which show Reynolds numbers on the order of 10 4 [START_REF] Warhaft | Turbulence in nature and in the laboratory[END_REF]. Hence DNS of turbulence will continue to be computationally demanding for the foreseeable future.

Due to the large number of grid points, practical implementations of DNS typically employ one-or two-dimensional domain decompositions within a distributed memory parallel programming paradigm. While the numerical solution of the field equations is typically achieved with well-established methods, the efficient implementation of particle tracking within such parallel approaches still poses major algorithmic challenges. In particular, particle tracking requires an accurate interpolation of the flow fields on distributed domains, and particles traversing the domain need to be passed on from one subdomain/process to another. As the Reynolds number increases, the number of particles required to adequately sample the turbulent fields needs to grow with the increasing number of grid points since this is a measure of the degrees of freedom of the flow. In addition higher-order statistics might be needed to address specific research questions, and thus the number of particles required for converged statistics increases as well [START_REF] Yeung | Lagrangian investigations of turbulence[END_REF][START_REF] Biferale | Multifractal statistics of Lagrangian velocity and acceleration in turbulence[END_REF][START_REF] Eyink | Stochastic flux freezing and magnetic dynamo[END_REF][START_REF] Eyink | Flux-freezing breakdown in high-conductivity magnetohydrodynamic turbulence[END_REF][START_REF] Biferale | Intermittency in the relative separations of tracers and of heavy particles in turbulent flows[END_REF][START_REF] Johnson | Large-deviation joint statistics of the finitetime Lyapunov spectrum in isotropic turbulence[END_REF][START_REF] Lalescu | Acceleration statistics of tracer particles in filtered turbulent fields[END_REF][START_REF] Lalescu | How tracer particles sample the complexity of turbulence[END_REF]. Overall, this requires an approach which handles the parallel implementation in an efficient manner for arbitrarily accurate methods. One option explored in the literature is the use of the high-level programming concept of coarrays, in practice shifting responsibility for some of the communication operations to the compiler, see [START_REF] Buaria | A highly scalable particle tracking algorithm using partitioned global address space (pgas) programming for extreme-scale turbulence simulations[END_REF]. The general solution that we describe makes use of MPI and OpenMP for explicit management of hardware resources. The combination of MPI [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message Passing Interface, 2nd Edition[END_REF] and OpenMP [START_REF] Board | OpenMP application program interface version 4.5[END_REF] has become a de facto standard in the development of large-scale applications [START_REF] Görler | The global version of the gyrokinetic turbulence code gene[END_REF][START_REF] Mininni | A hybrid MPI-OpenMP scheme for scalable parallel pseudospectral computations for fluid turbulence[END_REF][START_REF] Pekurovsky | P3DFFT: a framework for parallel computations of Fourier transforms in three dimensions[END_REF]42,[START_REF] Clay | A dual communicator and dual grid-resolution algorithm for petascale simulations of turbulent mixing at high schmidt number[END_REF][START_REF] Chatterjee | Scaling of a Fast Fourier Transform and a Pseudo-spectral Fluid Solver up to 196608 cores[END_REF]. MPI [START_REF] Walker | MPI: A standard message passing interface[END_REF] is used for communication between processes and OpenMP to manage multiple execution threads over multicore CPUs using shared memory. Separately, large datasets must be processed with specific data-access patterns to make optimal use of modern hardware, as explained for example in [START_REF] Homann | SoAx: A generic C++ Structure of Arrays for handling particles in HPC codes[END_REF].

To address the challenges outlined above, we have developed the numerical framework "Turbulence Tools: Lagrangian and Eulerian" (TurTLE), a flexible pseudo-spectral solver for fluid and turbulence problems implemented in C++ with a hybrid MPI/OpenMP approach [START_REF]TurTLE source code[END_REF][START_REF]TurTLE documentation[END_REF]. TurTLE allows for an efficient tracking of various types of particles. In particular, TurTLE showcases a parallel programming pattern for particle tracking that is easy to adapt and implement, and which allows for efficient use of available resources. Our eventdriven approach is especially suited for the case where individual processes require data exchanges with several other processes while also being responsible for local work. For this, asynchronous inter-process communication and tasks are used, based on a combined MPI/OpenMP implementation. As we will show in the following, TurTLE permits numerical particle tracking at relatively small costs, while retaining flexibility with respect to number of particles and numerical accuracy. We show that TurTLE scales well up to O(10 4 ) computing cores, with the flow field solver approximately retaining the performance of the used Fourier transform libraries for DNS with 3 × 2048 3 and 3 × 4096 3 degrees of freedom. We also measure the relative cost of tracking up to 2.2 × 10 9 particles as approximately only 10% of the total wall-time for the 4096 3 case, demonstrating the efficiency of the new algorithm even for very demanding particle-based studies.

In the following, we introduce TurTLE and particularly focus on the efficient implementation of particle tracking. Section 2 introduces the evolution equations for the fluid and particle models, as well as the corresponding numerical methods. Section 3 provides an overview of our implementation, including a more detailed presentation of the parallel programming pattern used for particle tracking. Finally, Section 4 summarizes a performance evaluation using up to 512 computational nodes.

Evolution equations and numerical method

Fluid equations

While TurTLE is developed as a general framework for a range of fluid equations, we focus on the Navier-Stokes equations as a prototypical example in the following. The incompressible Navier-Stokes equations take the form

∂ t u + u • ∇u = -∇p + ν∆u + f , ∇ • u = 0. (1)
Here, u(x, t) denotes the three-dimensional velocity field, p(x, t) is the kinematic pressure, ν is the kinematic viscosity, and F(x, t) denotes an external forcing that drives the flow. We consider periodic boundary conditions, which allows for the use of a Fourier pseudo-spectral scheme. Within this scheme, a finite Fourier representation is used for the fields, and the nonlinear term of the Navier-Stokes equations is computed in real space -an approach pioneered by Orszag and Patterson [START_REF] Orszag | Numerical simulation of three-dimensional homogeneous isotropic turbulence[END_REF].

For the concrete implementation in TurTLE, we use the vorticity formulation of the Navier-Stokes equation, which takes the form

∂ t ω = ∇ × (u × ω) + ν∆ω + F, (2)
where ω(x, t) = ∇ × u(x, t) is the vorticity field and F(x, t) = ∇ × f (x, t) denotes the curl of the Navier-Stokes forcing. The Fourier representation of this equation takes the form [START_REF] Homann | Lagrangian statistics of turbulent flows in fluids and plasmas[END_REF][START_REF] Wilczek | Statistical and numerical investigations of fluid turbulence[END_REF]

∂ t ω(k, t) = ik × F [u(x, t) × ω(x, t)] -νk 2 ω(k, t) + F(k, t), ( 3
)
where F is the direct Fourier transform operator. In Fourier space, the velocity can be conveniently computed from the vorticity using Biot-Savart's law,

û(k, t) = ik × ω(k, t) k 2 . ( 4
)
Equation ( 3) is integrated with a third-order Runge-Kutta method [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF], which is an explicit Runge-Kutta method with the Butcher tableau

0 1 1 1/2 1/4 1/4 1/6 1/6 2/3 . (5)
In addition to the stability properties described in [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF], this method has the advantage that it is memory-efficient, requiring only two additional field allocations, as can be seen from [START_REF] Homann | Lagrangian statistics of turbulent flows in fluids and plasmas[END_REF] ŵ1

(k) = ( ω(k, t) + hN[ ω(k, t)])e -νk 2 h , ŵ2 (k) = 3 4 ω(k, t)e -νk 2 h/2 + 1 4 ( ŵ1 (k) + hN[ ŵ1 (k)])e νk 2 h/2 , ω(k, t + h) = 1 3 ω(k, t)e -νk 2 h + 2 3 ( ŵ2 (k) + hN[ ŵ2 (k)])e -νk 2 h/2 , (6)
where h is the time step, limited in practice by the Courant-Friedrichs-Lewy (CFL) condition [START_REF] Courant | On the Partial Difference Equations of Mathematical Physics[END_REF]. The nonlinear term

N[ ŵ(k)] = ik × F F -1 ik× ŵ(k) k 2 × F -1 [ ŵ(k)] (7)
is computed by switching between Fourier space and real space.

If the forcing term is nonlinear, it can be included in the righthand side of [START_REF] Stohl | A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements[END_REF]. To treat the diffusion term, we use the standard integrating factor technique [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF] in [START_REF] Homann | Clustering of passive impurities in magnetohydrodynamic turbulence[END_REF].

Equation (3) contains the Fourier transform of a quadratic nonlinearity. Since numerical simulations are based on finite Fourier representations, the real-space product of the two fields will in general contain unresolved high-wavenumber harmonics, leading to aliasing effects [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. In TurTLE, de-aliasing is achieved through the use of a smooth Fourier filter, an approach that has been shown in [START_REF] Hou | Computing nearly singular solutions using pseudospectral methods[END_REF] to lead to good convergence to the true solution of a PDE, even though it does not completely re-move aliasing effects.

The Fourier transforms in TurTLE are evaluated using the FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. Within the implementation of the pseudospectral scheme, the fields have two equivalent representations: an array of vectorial Fourier mode amplitudes, or an array of vectorial field values on the real-space grid. For the simple case of 3D periodic cubic domains of size [0, 2π] 3 , the real-space grid is a rectangular grid of N × N × N points, equally spaced at distances of δ ≡ 2π/N. Exploiting the Hermitian symmetry of real fields, the Fourier-space grid consists of N × N × (N/2 + 1) modes. Therefore, the field data consists of arrays of floating point numbers, logically shaped as the real-space grid or arrays of floating point number pairs (e.g. fftw complex) logically shaped as the Fourier-space grid. Extensions to non-cubic domains or non-isotropic grids are straightforward.

The direct numerical simulation algorithm then has two fundamental constructions: loops traversing the fields, with an associated cost of O(N 3 ) operations, and direct/inverse Fourier transforms, with a cost of O(N 3 log N) operations.

Particle equations

A major feature of TurTLE is the capability to track different particle types, including Lagrangian tracer particles, ellipsoids, self-propelled particles and inertial particles. To illustrate the implementation, we focus on tracer particles in the following.

Lagrangian tracer particles are virtual markers of the flow field starting from the initial position x. Their position X evolves according to d dt X(x, t) = u(X(x, t), t), X(x, 0) = x.

(

) 8
The essential characteristic of such particle equations is that they require as input the values of various flow fields at arbitrary positions in space.

TurTLE combines multi-step Adams-Bashforth integration schemes (see, e.g., §6.7 in [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]) with spline interpolations [START_REF] Lalescu | Implementation of high order spline interpolations for tracking test particles in discretized fields[END_REF] in order to integrate the ODEs. Simple Lagrange interpolation schemes (see, e.g., §3.1 in [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]) are also implemented in Tur-TLE for testing purposes. There is ample literature on interpolation method accuracy, efficiency, and adequacy for particle tracking, e.g. [START_REF] Yeung | An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence[END_REF][START_REF] Lekien | Tricubic interpolation in three dimensions[END_REF][START_REF] Homann | Impact of the floating-point precision and interpolation scheme on the results of dns of turbulence by pseudospectral codes[END_REF][START_REF] Van Hinsberg | On the efficiency and accuracy of interpolation methods for spectral codes[END_REF][START_REF] Van Hinsberg | Optimal interpolation schemes for particle tracking in turbulence[END_REF]. The common feature of all interpolation schemes is that they can be represented as a weighted real-space-grid average of a field, with weights given by the particle's position. For all practical interpolation schemes, the weights are zero outside of a relatively small kernel of grid points surrounding the particle, i.e. the formulas are "local". For some spline interpolations, a non-local expression is used, but it can be rewritten as a local expression where the values on the grid are precomputed through a distinct global operation [START_REF] Yeung | An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence[END_REF] -this approach, for example, is used in [START_REF] Buaria | A highly scalable particle tracking algorithm using partitioned global address space (pgas) programming for extreme-scale turbulence simulations[END_REF].

Thus an interpolation algorithm can be summed up as follows:

1. compute X = X mod 2π (because the domain is periodic).

2. find the closest grid cell to the particle position X, indexed by c ≡ (c 1 , c 2 , c 3 ).

3. compute x = Xcδ.

4. compute a sum of the field over a kernel of I grid points in each of the 3 directions, weighted by some polynomials:

u(X) ≈ I/2 i 1 ,i 2 ,i 3 =1-I/2 β i 1 x1 δ β i 2 x2 δ β i 3 x3 δ u(c + i), (9)
where i = (i 1 , i 2 , i 3 ). The cost of the sum itself grows as I 3 , the cube of the size of the interpolation kernel. The polynomials β i j are determined by the interpolation scheme (see [START_REF] Lalescu | Implementation of high order spline interpolations for tracking test particles in discretized fields[END_REF]). Accuracy improves with increasing I. In TurTLE, interpolation is efficiently implemented even at large I. As discussed below in §3.3, this is achieved by sorting particle data such that a minimal number of MPI messages is used to complete all the triple sums, independently of the total number of particles.

Implementation

Overview

TurTLE follows the object-oriented paradigm to facilitate maintenance and extensions through class inheritance and virtualization. The solver relies on two types of objects. Firstly, all solvers inherit three basic elements from an abstract parent class: generic initialization, do work and finalization functionality. Solvers are then grouped in this first type of object, i.e. members of the resulting class hierarchy. The second type of object encapsulates essential data structures (i.e. fields, sets of particles) and associated functionality (e.g. HDF5-based I/O): these are "building block"-classes. Each solver then consists of a specific "arrangement" of the building blocks.

The parallelization of TurTLE is based on a standard, MPIbased, one-dimensional domain-decomposition approach: The three-dimensional fields are decomposed along one of the dimensions into a number of slabs, with each MPI process holding one such slab. In order to efficiently perform the costly FFT operations with the help of a high-performance numerical library such as FFTW, process-local, two-dimensional FFTs are interleaved with a global transposition of the data in order to perform the FFTs along the remaining dimension. A wellknown drawback of the slab decomposition strategy offered by FFTW is its limited parallel scalability, because at most N MPI processes can be used for N 3 data. We compensate for this by utilizing the hybrid MPI/OpenMP capability of FFTW (or functionally equivalent libraries such as Intel MKL), which allows to push the limits of scalability by at least an order of magnitude, corresponding to the number of cores of a modern multicore CPU or NUMA domain, respectively. All other relevant operations in the field solver can be straightforwardly parallelized with the help of OpenMP. Our newly developed parallel particle tracking algorithm has been implemented on top of this slab-type data decomposition using MPI and OpenMP, as shall be detailed below. Slab decompositions are beneficial for particle tracking since MPI communication overhead is minimized compared to, e.g., two-dimensional decompositions. Fields are split into slabs and distributed between P MPI processes along the x 3 direction. The N p particles are also distributed, with each MPI process storing P p particles on average. Within each MPI process the particle data is sorted according to its x 3 location. This leads to a direct association between each of the P s field slices to contiguous regions of the particle data arrays -in turn simplifying the interpolation procedure (see text for details). On average, S p particles are held within each such contiguous region.

Fluid solver

The fluid solver consists of operations with field data, which TurTLE distributes among a total of P MPI processes with a standard slab decomposition, see Fig. 1. Thus the logical field layouts consist of (N/P) × N × N points for the real-space representation, and (N/P)×N ×(N/2+1) points for the Fourier-space representation. This allows the use of FFTW [START_REF] Frigo | The design and implementation of FFTW3[END_REF] to perform costly FFT operations, as outlined above. We use the convention that fields are distributed along the real-space x 3 direction, and along the k 2 direction in the Fourier space representation (directions 2 and 3 are transposed between the two representations). Consequently, a problem on an N 3 grid can be parallelized on a maximum of N computational nodes using one MPI process per node and, possibly, OpenMP threads inside the nodes, see Fig. 1.

In the interest of simplifying code development, TurTLE uses functional programming for the costly traversal operation. Functional programming techniques allow to encapsulate field data in objects, while providing methods for traversing the data and computing specified arithmetic expressions -i.e. the class becomes a building block. While C++ allows for overloading arithmetic operators as a mechanism for generalizing them to arrays, our approach permits to combine several operations in a single data traversal, and it applies directly to operations between arrays of different shapes. In particular operations such as the combination of taking the curl and the Laplacian of a field (see [START_REF] Yeung | Lagrangian statistics from direct numerical simulations of isotropic turbulence[END_REF]) are in practice implemented as a single field traversal operation.

Particle tracking

We now turn to a major feature of TurTLE: the efficient tracking of particles. The novelty of our approach warrants a more in-depth presentation of the data structure and the parallel al-gorithms, for which we introduce the following notations (see also Fig. 1):

• P : the number of MPI processes (should be a divisor of the field grid size N);

• P s = N/P : the number of field slices in each slab;

• N p : the number of particles in the system;

• P p : the number of particles contained in a given slab (i.e. hosted by the corresponding process) -on average equal to N p /P;

• S p : the number of particles per slice, i.e. number of particles found between two slices -on average equal to N p /N;

• I : the width of the interpolation kernel, i.e. the number of slices needed to perform the interpolation.

The triple sum ( 9) is effectively split into I double sums over the x 1 and x 2 directions, the results of which then need to be distributed/gathered among the MPI processes such that the sum along the x 3 direction can be finalized. Independently of P and N, there will be N p sums of I 3 terms that have to be performed. However, the amount of information to exchange depends on the DNS parameters N, N p , and I, and on the job parameter P.

Whenever more than one MPI process is used, i.e. P > 1, we distinguish between two cases:

1. I ≤ P s , i.e. each MPI domain extends over at least as many slices as required for the interpolation kernel. In this case particles are shared between at most two MPI processes, therefore each process needs to exchange information with two other processes. In this case, the average number of shared particles is S p (I -1).

2. I > P s , i.e. the interpolation kernel always extends outside of the local MPI domain. The average number of shared particles is S p P s . Each given particle is shared among a maximum of I/P s processes, therefore each process must in principle communicate with 2 I/(2P s ) other processes.

The second scenario is the more relevant one for scaling studies.

Our expectation is that the communication costs will outweigh the computation costs, therefore the interpolation step should scale like N p I/P s = N p IP/N. In the worst case scenario, when the individual double sums have a significant cost as well, we expect scaling like N p I 3 P/N.

Particle data structure

The field grid is partitioned in one dimension over the processes, as described in Section 3.2, such that each process owns a field slab. For each process, we use two arrays to store the data for particles included inside the corresponding slab. The first array contains state information, including the particle locations -required to perform the interpolation of the field. We denote this first array state. The second array, called rhs, contains the value of the right-hand side of [START_REF] Stohl | Characteristics of atmospheric transport into the arctic troposphere[END_REF], as computed at the most recent few iterations (as required for the Adams-Bashforth integration); updating this second array requires interpolation. The two arrays use an array of structures pattern, in the sense that data associated to one particle is contiguous in memory. While this may lead to performance penalties, as pointed out in [START_REF] Homann | SoAx: A generic C++ Structure of Arrays for handling particles in HPC codes[END_REF], there are significant benefits for our MPI parallel approach, as explained below. We summarize in the following the main operations that are applied to the arrays.

Ordering the particles locally. When N > P, processes are in charge of more than one field slice, and the particles in the slab are distributed across different slices. In this case, we store the particles that belong to the same slice contiguously in the arrays, one slice after the other in increasing x 3 -axis order. This can be achieved by partitioning the arrays into P s different groups (indexed by s in the following) and can be implemented as an incomplete Quicksort with a complexity of O(P p log P s ) on average. The different groups will in general contain different numbers of particles, which implies that the associated data is stored at arbitrary locations within the contiguous state and rhs arrays. We therefore build an array offset of size P s + 1, where offset[s] returns the starting index of the first particle for the group s and the number of particles in this same group s can be computed as offset[s+1]-offset [s]. For simplicity, we explicitly store offset[P s ]= P p as the P s + 1-th element of offset. This allows direct access to the contiguous data regions corresponding to each field slice, access which is in turn relevant for MPI exchanges (see below).

Exchanging the particles for computation. With our data structures, we are able to send the state information of all the particles located in a single group with only one communication, which reduces communication overhead. Moreover, sending the particles from several contiguous groups can also be done in a single operation because the groups are stored sequentially inside the arrays.

Particles displacement/update. The positions of the particles are updated at the end of each iteration, and so the arrays must be rearranged accordingly. The changes in the x 3 direction might move some particles to a different slice and even to a slice owned by a different process. Therefore, we first partition the first and last groups (the groups of the first and last slices of the process's slab) to move the particles that are now outside of the process's grid interval at the extremities of the arrays. We only act on the particles located at the lower and higher groups because we assume that the particles cannot move a distance greater than 2π/N. For regular tracers [START_REF] Stohl | Characteristics of atmospheric transport into the arctic troposphere[END_REF] this is in fact required by the CFL stability condition of the fluid solver. This partitioning is done with a complexity O(P p /P s ). Then, every process exchanges those particles with its direct neighbors, ensuring that the particles are correctly distributed. Finally, each process sorts its particles to take into account the changes in the positions and the newly received particles as described previously.

Parallelization

The interpolation of the field at the particle locations concentrates most of the workload of the numerical particle tracking. For each particle, the interpolation uses the I 3 surrounding field nodes. However, because we do not mirror the particle or the field information on multiple processes, we must actively exchange either field or particle information to perform a complete interpolation. Assuming that the number of particles in the simulation is much smaller than the number of field nodes, i.e. the relation P p < IN 2 holds, less data needs to be transferred on average when particle locations are exchanged rather than field values at required grid nodes. Consequently, in our implementation we exchange the particle information only.

A straightforward implementation in which the communication and computation are dissociated consists in the following operations:

(a) each process computes the interpolation of its particles on its field;

(b) all the processes exchange particle positions with their neighbors (each process sends and receives arrays of positions);

(c) each process computes the interpolation using its field on the particle positions it received from other processes in (b);

(d) all the processes exchange the results of the interpolations from (c) with the corresponding neighbors;

(e) each process merges the results it received in (d) and the results from its own computation from (a).

In our implementation, we interleave these five operations to overlap communication with computation. As we detail in the following, the master thread of each MPI process creates computation work packages, then performs communications while the other threads are already busy with the work packages. This is achieved with the use of non-blocking MPI communications and OpenMP tasks, as illustrated in Fig. 2. Broadly speaking, the execution is guided by a pairing of communications with relevant actions that need to be performed on the associated data (either sent or received), with all such pairs put together in a list L (cap list in Fig. 2). Actions of "compute" type are turned into tasks distributed to all available threads, whereas actions of "send/receive" type are handled separately by the master thread. A more detailed overview follows. The figure shows a simplified version of the parallel execution pattern (the function exchange compute), along with a possible execution timeline. For simplicity we also define the function proceed, which encompasses much of the communication and coordination activity of the master thread. Here we focus on describing the algorithm, only referring to the execution timeline implicitly through the color correspondence. In a first stage, the master thread splits the local interpolation from (a) into tasks and submits them immediately but with a low priority. Then, it posts all the sends/receives related to (b) and all the receives related to (d), and creates the list L. In the core part of the algorithm, the master thread performs a wait-any on the list of MPI requests (line 7 in the pseudocode). This MPI function is blocking and returns as soon as one of the communications in the list is completed. Hence, when a communication is completed, the master thread acts accordingly to the type of action A it represents (a.type in Fig. 2). If A is the completion of a send of local particle positions, from (b), there is nothing to do and the master thread directly goes back to the wait-any. In this case, it means that a send is completed and that there is nothing new to do locally. If A is the completion of a receive of remote particle positions, from (b), then the master thread creates tasks to perform the interpolation of these positions, from (c), and submits them with high priority (lines 23 and 24 in the pseudocode). Setting a high priority ensures that all the threads will work on these tasks even if the tasks inserted earlier to interpolate the local positions, from (a), are not completed. When these tasks are completed, the master thread posts a non-blocking send to communicate the results to the process that owns the particles and adds the corresponding MPI request to L. Then, the master thread goes back to the wait-any on L. If A is the completion of a send of interpolation on received positions, as just described, the master thread has nothing to do and goes back to the waitany on L. In fact, this event simply means that the results were correctly sent. If A is the completion of a receive, from (d), of interpolation performed by another process, done in (c), the master thread keeps the buffer for merging at the end, and it goes back to the wait-any. When L is empty, it means that all communications (b,d) but also computations on other positions (c) are done. If some local work still remains from (a), the master thread can join it and compute some tasks. Finally, when all computation and communication are over, the threads can merge the interpolation results, operation (e) (line 17 in the pseudocode). The described strategy is a parallel programming pattern that could be applied in many other contexts when there are local and remote works to perform and where remote work means first to exchange information and second to apply computation on it.

In-Order Parallel Particle I/O

Saving the states of the particles on disk is a crucial operation to support checkpoint/restart and for post-processing. We focus on output because TurTLE typically performs many more output than input operations (the latter only happen during initialization). The order in which the particles are saved is important because it influences the writing pattern and the data accesses during later post-processing of the files. As the particles move across the processes during the simulation, a naive output of the particles as they are distributed will lead to inconsistency from one output to the other. Such a structure would require reordering the particles during the post-processing or would result in complex file accesses. That is why we save the particles in order, i.e. in the original order given as input to the application.

The algorithm that we use to perform the write operation is shown in Fig. 3. There are four main steps to the procedure: pre-sort ("Sort" and "Find pivot" in the figure), followed by exchange ("Send/Recv" in Fig. 3), with a final post-sort ("Sort") before the actual HDF5 write ("Parallel write to file" in the figure).

Each process first sorts its local particles using the global indices with O(P p log P p ) complexity. This sort can be done in Here, we consider that the particle data (illustrated by the global particle index) is distributed among 4 processes, but that only 2 of them are used in the write operation. parallel using multiple threads. Then, each process counts the number of particles it has to send to each of the processes that are involved in the file writing. These numbers are exchanged between the processes allowing each process to allocate the reception buffer. If we consider that P O processes are involved in the output operation, each of them should receive N p /P O particles in total, and a process of rank r should receive the particles from index r × N p /P O to (r + 1) × N p /P O -1.

In the exchange step, the particles can be sent either with multiple non-blocking send/receive operations or with a single all-to-all operation, with the total number of communications bounded by P × P O . Finally, the received particles are sorted with a complexity of O(N p /P O log N p /P O ), and written in order into the output file.

The number P O of processes involved in the writing should be carefully chosen because as P O increases, the amount of data output per process decreases and might become so small that the write operation becomes inefficient. At the same time, the preceding exchange stage becomes more and more similar to a complete all-to-all communication with N 2 p relatively small messages. On the other hand, as P O decreases, the size of the messages exchanged will increase, and the write operation can eventually become too expensive for only a few processes, which could also run out of memory. This is why we heuristically fix P O using three parameters: the minimum amount of data a process should write, the maximum number of processes involved in the write operation, and a chunk size. As N p increases, P O increases up to the given maximum. If N p is large enough, the code simply ensures that P O -1 processes output the same amount of data (being a multiple of the chunk size), and the last process writes the remaining data. In our implementation, the parameters are chosen empirically (based on our experience with several HPC clusters running the IBM GPF-S/SpectrumScale parallel file system), and they can be tuned for specific hardware configurations if necessary.

We use a similar procedure for reading the particle state: P O processes read the data, they sort it according to spatial location, then they redistribute it to all MPI processes accordingly.

Computational performance

Hardware and software environment

To evaluate the computational performance of our approach, we perform benchmark simulations on the HPC cluster SuperMUC-NG from the Leibniz Supercomputing Centre (LRZ): we use up to 512 nodes containing two Intel Xeon Platinum 8174 (Skylake) processors with 24 cores each and a base clock frequency of 3.1 GHz, providing 96 GB of main memory. The network interconnect is an Intel OmniPath (100 Gbit/s) with a pruned fat-tree topology that enables non-blocking communications within islands of up to 788 nodes. We use the Intel compiler 19.0, Intel MPI 2019.4, HDF5 1.8.21 and FFTW 3.3.8. For our benchmarks, we always fully populate the nodes, i.e. the combination of MPI processes per node (ppn) and OpenMP threads per MPI process is chosen such that their product equals 48, and that the threads spawned by an MPI rank are confined within the NUMA domain defined by a single processor. 1: Execution times for TurTLE on SuperMUC-NG (in seconds), various configurations (values correspond to Fig. 4a).

Overall performance

Figure 4 provides an overview of the overall parallel scaling behavior of the code for a few typical large-scale setups (panel a as well as Table 1) together with a breakdown into the main algorithmic steps (panel b). We use the execution time for a single time step (averaged over a few steps) as the primary performance metric and all data and computations are handled in double precision. The left panel shows, for two different setups (N = 2048, 4096), that the code exhibits excellent strongscaling efficiency (the dashed line represents ideal scaling) from the minimum number of nodes required to fit the code into memory up to the upper limit which is given by the maximum number of MPI processes that can be utilized with our onedimensional domain decomposition. The total runtime of Tur-TLE, and hence its major parallel scaling properties are largely determined by the three-dimensional fast Fourier transforms shows a breakdown of the total runtime into the main functional elements, namely solving the system of Navier-Stokes partial differential equations ("PDE misc" and "PDE FFT") and particle tracking ("IFT for PT" and "PT"). The fluid solver cost is largely dominated by the fast Fourier transforms ("PDE FFT"). The cost of particle tracking for 10 8 particles (with I = 8) is determined by the additional inverse Fourier transform ("IFT for PT"), whereas the runtime for our novel particle tracking algorithm ("PT") is still negligible for 10 8 particles. Hatched regions represent the fraction of MPI communication times.

(using FFTW) for solving the system of Navier Stokes partial differential equations, as shown in Fig. 4b (labeled "PDE FFT", entire blue area) for the example of the large setup (N = 4096) with 8 processes per node (see also Table 2). With increasing node count, the runtime of the FFTs in turn gets increasingly dominated by an all-to-all type of MPI communication pattern which is arising from the FFTW-internal global transpositions (blue-hatched area) of the slab-decomposed data. While a good OpenMP efficiency can be noted for the case of 64 nodes by comparing the blue and the orange curve in Fig. 4a,i.e. the same problem computed with a different combination of MPI processes per node [START_REF] Stohl | Characteristics of atmospheric transport into the arctic troposphere[END_REF][START_REF] Durham | Turbulence drives microscale patches of motile phytoplankton[END_REF] and a corresponding number of OpenMP threads per process (6, 3), we generally observe FFTW exhibiting limited OpenMP scaling beyond a number of 6 to 8 threads per MPI process for the grid resolutions (N) considered here. Moreover, even at moderate thread counts, FFTW, due to internal load-imbalances, fails to efficiently handle the process-local transforms for certain dimensions of the local data slabs: In the case of 256 nodes, FFTW apparently cannot efficiently use more than 3 OpenMP threads for parallelizing over the local slabs of dimension 2 × 4096 × 2049 (Fourier representation), whereas good scaling up to the desired maximum of 6 threads is observed for a dimension of 8 × 4096 × 2049 (128 nodes) and also 1 × 4096 × 2049 (512 nodes). The same argument applies for the smaller setup (N = 2048) on 128 nodes and is the root cause for the kink in the strong scaling curves of panel a. We plan for TurTLE to support FFTs also from the Intel Math Kernel Library (MKL) which are expected to deliver improved threading efficiency. For practical applications, this is less of a concern, since a user can perform a few exploratory benchmarks for a given setup of the DNS on the particular computational platform, and available node counts, in order to find an optimal combination of MPI processes and OpenMP threads. Since the runtime per time step is constant for our implementation of the Navier-Stokes solver, a few time steps are sufficient for tuning a long-running DNS.

Thanks to our efficient and highly scalable implementation of the particle tracking, its contribution to the total runtime is barely noticeable in the figure ("PT", red colour in Fig. 4b). This holds even for significantly larger numbers of particles than the value of N p = 10 8 which was used here (see below for an in-depth analysis). The only noticeable additional cost for particle tracking, amounting to roughly 10% of the total runtime, comes from an additional inverse FFT ("IFT for PT", green colour) which is required to compute the advecting velocity field, which is independent of N p and scales very well.

Finally, Fig. 4a also suggests good weak scaling behavior of TurTLE: When increasing the resolution from N = 2048 to N = 4096 and at the same time increasing the number of nodes from 64 to 512, the runtime increases from 10.35s to 11.45s, which is consistent with a O(N 3 log N) scaling of the runtime, given the dominance of the FFTs.

Particle tracking performance

Fig. 5 provides an overview and some details of the performance of our particle tracking algorithm, extending the assess- ment of the previous subsection to particle numbers beyond the current state of the art [START_REF] Buaria | A highly scalable particle tracking algorithm using partitioned global address space (pgas) programming for extreme-scale turbulence simulations[END_REF]. We use the same setup of a DNS with N = 4096 and 8 MPI processes per node on SuperMUC-NG, as presented in the previous subsection. Fig. 5a summarizes the strong-scaling behavior on 128, 256 or 512 nodes of SuperMUC-NG for different numbers of particles (10 6 , 10 8 and 2.2 × 10 9 ) and for different sizes of the interpolation kernel I (4, 6, 8). Most importantly, the absolute run times are small compared to the fluid solver: Even when using the most accurate interpolation kernel, a number of 2.2 × 10 9 particles can be handled within less than a second (per time step), i.e. less than 10% to the total computational cost of Tur-TLE on 512 nodes per time step (cf. Fig. 4).

The case of N p = 10 6 is shown only for reference here. This number of particles is too small to expect good scalability in the regime of 128 compute nodes and more. Still, the absolute runtimes are minuscule compared to a DNS of typical size. For N p = 10 8 we observe good but not perfect strong scaling, in particular for the largest interpolation kernel (I = 8), suggesting that we observe the N p IP/N regime, as discussed previously. It is worth mentioning that we observe a sub-linear scaling of the total runtime with the total number of particles (Fig. 5a).

Fig. 5b shows a breakdown of the total runtime of the particle tracking algorithm into its main parts, interpolation (operations detailed in Fig. 2, shown in orange) and redistribution (local sorting of particles together with particle exchange, blue), together with the percentage of time spent in MPI calls. The latter takes between half and two thirds of the total runtime for N p = 10 8 particles (cf. upper panel b) and reaches almost 80% for N p = 2.2×10 9 particles on 512 nodes (lower panel c). Overall, the interpolation cost always dominates over redistribution, and increases with the size of the interpolation kernel roughly as I 2 , i.e. the interpolation cost is proportional to the number of MPI messages required by the algorithm (as detailed above). Figure 6 provides an overview of the computational costs of the main parts of the output algorithm, namely sorting particles according to their initial order (pre-sort and post-sort stages, cf. Sect. 3.4), communicating the data between processes (exchange stage), and writing data to disk using parallel HDF5 (write stage). Here, the same setup is used as in Fig. 5 panels b and c, respectively, noting that the output algorithm does not depend on the size of the interpolation kernel. The figure shows that the total time is largely dominated by the write and exchange stages, with the sorting stages not being significant. Of the latter, the post-sort operation is relatively more expensive than the pre-sort stage, because only a comparably small subset of P O < P processes is used in the post-sort stage (in the present setup P O = 1 for 10 6 particles, P O = 72 for 10 8 particles, and P O = 126 for 2.2 × 10 9 particles were used). This indicates that our strategy of dumping the particle data in order adds only a small overhead, which is mostly spent in the communication stage (unsorted output could be done with a more simple communication pattern) but not for the actual (processlocal) reordering of the particles. For a given number of particles N p , the number of processes P O involved in the write operation is fixed, independent of the total number P of processes used for the simulation. Consequently, the time spent in the write stage does not depend on the number of nodes (and hence P), as shown in Fig. 6a. However, P O may increase with increasing N p (and fixed P).

Fig. 6b shows that the cost of writing 10 6 particles with a single process is negligible, whereas writing 10 8 particles with 72 processes becomes significant, even though a similar number of particles per output process (1.4 × 10 6 particles) is used. This reflects the influence of the (synchronization) overhead of the parallel HDF5 layer and the underlying parallel IO system. On the other hand, it takes about the same amount of time for 126 processes to write 1.7 × 10 7 particles each, compared with 72 processes writing 1.4 × 10 6 particles each, which motivates our strategy of controlling the number of processes P O that are involved in the interaction with the IO system. However, the choice of P O also influences the communication time spent in the exchange stage. When looking at the exchange stage in Figure 6a, we recall that 72 processes write the data for all three node counts. As P increases, the 72 processes receive less data per message but communicate with more processes. From these results it appears that this is beneficial: reducing the size of the messages but increasing the number of processes that communicate reduces the overall duration of this operation (that we do not control explicitly since we rely on the MPI Alltoallv collective routine). For a fixed number of processes and an increasing number of particles (see Figure 6b), the total amount of data exchanged increases and the size of the messages varies. The number P O (i.e., 1, 72 and 126) is not increased proportionally with the number of particles N p (i.e., 10 6 , 10 8 and 2.2 × 10 9 ), which means that the messages get larger and, more importantly, each process needs to send data to more output processes. Therefore, increasing P O also increases the cost of the exchange stage but allows to control the cost of write stage. Specifically, it takes about 4s to output 10 8 particles (1s for exchange and 3s for write). It takes only 6 times longer, about 23s (15s for exchange, 4s for write, and 3s post-sort) to output 22 times more, 2.2 × 10 9 , particles.

Overall, our strategy of choosing the number of processes P O participating in the IO operations independent of the total number P of processes allows us to avoid performance-critical situations where too many processes would access the IO system, or too many processes would write small pieces of data. The coefficients used to set P O can be adapted to the specific properties (hardware and software stack) of an HPC system.

Summary and conclusions

In the context of numerical studies of turbulence, we have presented a novel particle tracking algorithm using an MPI/OpenMP hybrid programming paradigm. The implementation is part of TurTLE, which uses a standard pseudo-spectral approach for the direct numerical simulation of turbulence in a 3D periodic domain. TurTLE succeeds at tracking billions of particles with a negligible cost relative to solving the fluid equations. MPI communications are overlapped with computation thanks to a parallel programming pattern that mixes OpenMP tasks and MPI non-blocking communications. At the same time, the use of a contiguous and slice-ordered particle data storage allows to minimize the number of required MPI messages for any size of the interpolation kernel. This way, our approach combines both numerical accuracy and computational performance to address open questions regarding particle-laden flows by performing highly resolved numerical simulations on large supercomputers [START_REF] Lalescu | How tracer particles sample the complexity of turbulence[END_REF][START_REF] Bentkamp | Persistent accelerations disentangle lagrangian turbulence[END_REF][START_REF] Pujara | Shape-and scale-dependent coupling between spheroids and velocity gradients in turbulence[END_REF][START_REF] Bentkamp | The statistical geometry of material loops in turbulence[END_REF]. Indeed, as the benchmarks presented in Section 4 have shown, TurTLE scales up to many thousands of CPU cores on modern high-performance computers.

We expect that due to our task-based parallelization and the asynchronous communication scheme the particle-tracking algorithm is also well suited for offloading to massively parallel accelerators (e.g. GPUs) of a heterogeneous HPC node architecture. Concerning the fluid solver, the pseudo-spectral approach taken by TurTLE appears in principle amenable to GPU acceleration, as demonstrated by similar fluid codes (e.g. [START_REF] Mukherjee | Three dimensional pseudo-spectral compressible magnetohydrodynamic GPU code for astrophysical plasma simulation[END_REF][START_REF] Ravikumar | GPU acceleration of extreme scale pseudo-spectral simulations of turbulence using asynchronism[END_REF][START_REF] López | nsCouette -A high-performance code for direct numerical simulations of turbulent Taylor-Couette flow[END_REF][START_REF] Rosenberg | GPU parallelization of a hybrid pseudospectral geophysical turbulence framework using CUDA[END_REF]). Whether it can be implemented in a (performance-) portable and scalable way on large GPU-accelerated HPC systems remains to be investigated.

Figure 1 :

 1 Figure 1: Distribution of real-space data between MPI processes in TurTLE.Fields are split into slabs and distributed between P MPI processes along the x 3 direction. The N p particles are also distributed, with each MPI process storing P p particles on average. Within each MPI process the particle data is sorted according to its x 3 location. This leads to a direct association between each of the P s field slices to contiguous regions of the particle data arrays -in turn simplifying the interpolation procedure (see text for details). On average, S p particles are held within each such contiguous region.

1 function exchange_compute() 2 #pragma omp parallel 3 #pragma omp master 4 /

 234 / Create list of communication-action pairs (cap_list).

5 cap_list ← create_first_coms() 6 while cap_list is not empty do 7 index_done 1 Figure 2 :

 56712 Figure2: Overview of parallel interpolation algorithm (left) and reference execution timeline for one process with four threads (right). Colored lines of code correspond to individual timeline elements. The threads compute local operations by default but switch to remote operations when the master thread creates the corresponding new higher-priority tasks. With the use of priorities, the execution pattern allows for quicker communication of results to different processes. See text for details.

Figure 3 :

 3 Figure3: The different stages to perform a parallel saving of the particles in order. Here, we consider that the particle data (illustrated by the global particle index) is distributed among 4 processes, but that only 2 of them are used in the write operation.

Figure 4 :

 4 Figure 4: Computational performance of TurTLE. Strong scaling behavior of the total runtime (a) for grids with N = 2048 and N = 4096, respectively, and using 8 or 16 MPI processes per node (ppn) on up to 512 fully populated nodes (24576 cores) of SuperMUC-NG (I = 8, N p = 10 8 ). For N = 4096 and 8 ppn panel (b)shows a breakdown of the total runtime into the main functional elements, namely solving the system of Navier-Stokes partial differential equations ("PDE misc" and "PDE FFT") and particle tracking ("IFT for PT" and "PT"). The fluid solver cost is largely dominated by the fast Fourier transforms ("PDE FFT"). The cost of particle tracking for 10 8 particles (with I = 8) is determined by the additional inverse Fourier transform ("IFT for PT"), whereas the runtime for our novel particle tracking algorithm ("PT") is still negligible for 10 8 particles. Hatched regions represent the fraction of MPI communication times.

a

  PT execution times on SuperMUC-NG I = 4, N p = 10 6 I = 6, N p = 10 6 I = 8, N p = 10 6 I = 4, N p = 10 8 I = 6, N p = 10 8 I = 8, N p = 10 8 I = 4, N p = 2.2×10 9 I = 6, N p = 2.2×10 9 I = 8, N p = 2.2×10

8 I = 4 I = 6 I = 8 I = 4 I = 6 I = 8 I = 4 I = 6 Figure 5 :

 84848465 Figure 5: Computational performance of the particle tracking code using 8 MPI processes per node and a DNS of size N = 4096, for different sizes of the interpolation kernel I. Panel (a): strong scaling for different numbers of particles N p and sizes of the interpolation kernel (memory requirements limit the N p = 2.2 × 10 9 case to 512 nodes). The dashed line corresponds to ideal strong scaling. Panel (b): contributions of interpolation and redistribution operations to the total execution time, for a fixed number of particles, N p = 10 8 , and for different sizes of the interpolation kernel (the corresponding vertical bars are distinguished by hatch style, see labels on top) as a function of the number of compute nodes. Panel (c): relative contributions of interpolation and redistribution as a function of N p . Percentages represent the fraction of time spent in MPI calls.

Figure 6 :

 6 Figure6: Performance of particle output, as distributed between the four different operations (see also Fig.3): pre-sort (particles are sorted by each process), MPI exchange (particle data is transferred to processes actually participating in I/O), post-sort (particles are sorted on each I/O process), and write (HDF5 write call). Panel (a): elapsed times as a function of the total number of nodes, for a fixed N p = 10 8 (see also Fig.5b). Panel (b): elapsed time as a function of the number of particles, in the case of 512 nodes (see also Fig.5c).

Table

		N = 2048	N = 4096
	ppn
	Number	8	16	8	16
	of nodes
	64	10.3 9.43	-	-
	128	6.92 4.76 41.6	-
	256	2.52	-	31.1 17.5
	512	-	-	11.9	-

Table 2 :

 2 Breakdown of the total TurTLE runtime (in seconds) into main functional elements (values correspond to Fig.4b).

Acknowledgments

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (www.lrz.de). Some computations were also performed at the Max Planck Computing and Data Facility. This work was supported by the Max Planck Society. This work is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 101001081).


Jonatha Anselmi
email: jonatha.anselmi@inria.fr

Bruno Gaujal
email: bruno.gaujal@inria.fr

Energy Optimal Activation of Processors for the Execution of a Single Task with Unknown Size

Keywords: Energy minimization, GPU cluster, optimal speed, convex programming, parallel computing

A key objective in the management of modern computer systems consists in minimizing the electrical energy consumed by processing resources while satisfying certain target performance criteria. In this paper, we consider the execution of a single task with unknown size on top of a service system that offers a limited number of processing speeds, say N , and investigate the problem of finding a speed profile that minimizes the resulting energy consumption subject to a deadline constraint. Existing works mainly investigated this problem when speed profiles are continuous functions. In contrast, the novelty of our work is to consider discontinuous speed profiles, i.e., a case that arises naturally when the underlying computational platform offers a finite number of speeds.

In our main result, we show that the computation of an optimal speed profile boils down to solving a convex optimization problem. Under mild assumptions, for such convex optimization we prove some structural results that yield the formulation of an extremely efficient solution algorithm. Specifically, we show that the optimal speed profile can be computed by solving O(log N ) one-dimensional equations. Our results hold when the task size follows a known probability distribution function and the set of available speeds, if listed in increasing order, forms a sublinear concave sequence.

I. INTRODUCTION

To handle the demand induced by current deep learning technologies and applications, which affect several aspects of our daily life, AI companies and cloud providers are increasingly relying on large-scale GPU clusters. The electrical energy consumed by such systems has become a dominant issue because of its enormous financial and environmental impact [START_REF] Beloglazov | Optimal online deterministic algorithms and adaptive heuristics for energy and performance efficient dynamic consolidation of virtual machines in cloud data centers[END_REF]- [START_REF] Borgetto | Energy-aware service allocation[END_REF]. In this context, a key objective consists in minimizing the energy consumed by processing resources such as CPUs or GPUs while satisfying some target user-perceived performance criteria. In this paper, we consider the execution of a single task on top of a processing system composed of N processors and investigate how to design a scheduler that dynamically activates the optimal number of processors at any point in time. Here, optimality refers to minimizing the overall electrical energy consumed by executing the task subject to a deadline constraint. In cloud and super computing, distributing a task over multiple processors is a common practice. Within This work is supported by the French National Research Agency in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02) and the LabEx PERSYVAL (ANR-11-LABX-0025-01). machine learning tasks like TensorFlow, which are highly parallel, the number of processors used by a single task can vary by five orders of magnitude. This has been recently confirmed in a real trace from Google's Borg scheduler [START_REF] Tirmazi | Borg: The next generation[END_REF]. We also consider the realistic case where i) the processor speed-up function is sublinear and concave [START_REF] Zhan | Parsec3.0: A multicore benchmark suite with network stacks and splash-2x[END_REF], and ii) the scheduler does not know in advance the exact size of the task, though we assume that it knows its probability distribution function. The N -processors assumption is a mere modeling abstraction as our framework only requires that the underlying computational platform offers a finite number of speeds. This case can be met even within a single CPU provided that it is compatible with Dynamic Voltage Frequency Scaling (DVFS) technologies.

A. Relevant Related Work

In the literature, the problem of minimizing the energy consumption induced by task executions while satisfying realtime constraints has been investigated in several contexts depending on the information available to the scheduler and the type of computing resources. Without being exhaustive, our goal here is to provide a high-level review of the state of the art highlighting the difference with respect to our work.

The work closest to our own is [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF], which considers the execution of a single task as we do but when i) the speeds available to the scheduler form an interval of the real numbers and ii) the power consumption under speed s is of the form P (s) = s α with α > 2. Within these assumptions, the optimal speed profile admits an elegant "closed-form" expression. The "continuous speeds" assumption is justified by DVFS technologies as these allow the scheduler to change the current speed by changing the current CPU voltage. In practice however, processors may offer only a fixed, limited set of valid speeds, and in this case the solution in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF] does not extend to discrete speeds (as also mentioned in that reference). However, the following sub-optimal heuristic is proposed to handle in this case: "solve the continuous case and then use the closest discrete speed available to the scheduler". Unfortunately, it is not known how to bound the gap to the optimal solution, as also specified in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF]Section 8.4]. An additional difference of our work with respect to [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF] is that we put milder assumptions on the structure of the power consumption as a function of the processing speed; our framework includes functions of the form P (s) = s α but also all increasing convex functions such that P (0) = 0.

Other relevant related scenarios consider the case where the given task is composed of a stream of sporadic jobs all with known sizes, arrival times and deadlines. Here, the problem of designing an algorithm for the optimal speed profile, or the number of active processors to use at any point in time, has been solved in a series of papers, with improving complexity [START_REF] Gaujal | A Pseudo-Linear Time Algorithm for the Optimal Discrete Speed Minimizing Energy Consumption[END_REF], [START_REF] Li | An O(n 2 ) algorithm for computing optimal continuous voltage schedules[END_REF], [START_REF] Yao | A scheduling model for reduced CPU energy[END_REF]. In this online clairvoyant case, where the size and deadline of each job are revealed at the arrival time of the task, the problem can be formulated as a Markov Decision process and has been solved in [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF] resp. [START_REF] Anselmi | Optimal Speed Profile of a DVFS Processor under Soft Deadlines[END_REF] when the deadlines are hard resp. soft.

B. Contribution

Within the framework described above, our main contribution is to show that optimal speed profiles, or equivalently the number of processors to use at any point in time, are given by the solution of a convex optimization problem (see Theorem 1). This holds true when the task size follows an arbitrary probability distribution function and when the processor speed-up function is sublinear and concave. These speed-up function assumptions are common and justified because the process of distributing a single task over multiple processors usually implies some resource contention phenomena, and this slows down the computation [START_REF] Zhan | Parsec3.0: A multicore benchmark suite with network stacks and splash-2x[END_REF]. When the distribution of the task size is continuous and power consumption function is strictly convex, we use this convex optimization problem to show that optimal speed profiles possess the following structural properties: i) there exists a unique optimal speed profile, ii) the optimal speed profile is increasing, which means that the task execution accelerates over time, and iii) the optimal speed profile uses a set of consecutive speeds. Property ii) is intuitive and also holds within the continuousspeed framework considered in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF]. Property iii) yields the formulation of an extremely efficient solution algorithm (see Algorithm 1). Specifically, we show that the optimal speed profile can be computed by solving O(log N ) one-dimensional equations. The proposed algorithm is particularly useful to speed up the computation of an optimal profile in the context of machine learning tasks, as the number of processors used in this setting, and thus the number of available speeds, can vary by five orders of magnitude [START_REF] Tirmazi | Borg: The next generation[END_REF], e.g., N ∈ {1, . . . , 10 5 }.

C. Organization

This paper is organized as follows. In Section II, we describe our framework and define our energy minimization problem. In Section III, we show that the considered energy minimization problem belongs to the framework of convex programming. In Section IV, we propose Algorithm 1 for the efficient computation of the optimal speed profile. Finally, in Section VI, we draw the conclusions of our work and discuss some generalizations of our results.

II. ENERGY MINIMIZATION FRAMEWORK

The energy minimization framework that we investigate in the following is composed of the following components.

A. Processors

The underlying computational platform is composed of a finite number of identical processors, say N . Each processor operates with unitary speed, e.g., one work unit per second. By scaling work units, this assumption is not a loss of generality. When active resp. inactive, each processor consumes power P on , resp. P off . In other words, if a processor is active in a time interval [t, t ], then its energy consumption is (t -t)P on .

B. Tasks

A task is a sequence of unitary operations of random size W (the total amount of work) that arrives at time zero. Only statistical information are known about the task size W . In particular, it is known its size probability distribution function F (w) := P(W ≤ w), the probability that the task size is no more than w. Let also F c (w) := P(W > w). We allow F to be discontinuous, though in some of our results we will require F to be continuous, which implies that W admits a density function. The continuity assumption will simplify the numerical computation of an optimal speed profile. We also assume that the support of W is [W min , W max ], where W min ≥ 0 and W max < ∞ are the minimum and maximum task sizes, respectively. Finally, we impose the following hard real-time constraint: the task must be completed before a deadline D. This deadline must satisfy a feasibility constraint that will be discussed below.

C. Scheduler

The given task can be executed on multiple processors at the same time. The role of the scheduler is to allocate processors to the execution of the task at any point in time. The number of active processors may be changed at any time and we assume that these changes are immediate and do not incur any additional energy cost; for instance, ramp-up and/or speed change effects are neglected. The information available to the scheduler is:

• The deadline D (the task must be completed before D); • The task size distribution (with support [W min , W max ]);

• The effective speed s n of execution when n processors are active, for all n = 1, . . . , N ; • The power dissipation P (s n ) when n processors are active, for all n = 1, . . . , N ; • The current executed work w. The whole process stops as soon as the task is complete and all processors are put to rest immediately. Of course, the scheduler does not know when this is going to happen since it does not know the actual size of the task W .

We denote by s n = dw/dt the overall processing speed when allocating n processors to the task. The function s n is often called the speed-up function of the task and measures its scalability, or degree of parallelism. We assume that (s n ) n is non-decreasing and that s n+1 -s n is non-increasing. Thus, at higher speed, the work performed by unit of energy is lower, i.e., more energy is required to perform a fixed amount of work. These "diminishing returns" assumption is natural to model the overheads induced by parallel computations [START_REF] Zhan | Parsec3.0: A multicore benchmark suite with network stacks and splash-2x[END_REF].

We assume that the (instantaneous) power consumption when operating at speed s n , for n = 1, . . . , N , is

P (s n ) := nP on +(N -n)P idle = n(P on -P idle )+N P idle , (1)
where P on is the (instantaneous) power consumption of a processor turned on and P idle is the power consumption when the node is idle. We let P idle ≤ P on , and thus P (s n ) is increasing in n. We also assume P (0) = 0, which is equivalent to assume that the system is at complete rest when no processor is active.

Let also P ∆ = P on -P idle and P off = N P idle . Given that the speed-up sequence (s n ) n has non-increasing increments, P (s) is convex as long as s 2 (P ∆ + P off ) ≤ s 1 (2P ∆ + P off ), with strict inequality holding as long as P (s) is strictly convex. The latter inequality, which will be implicitly assumed in the following, holds true as long as P off is not too large and/or the speed-up is low enough.

Since the scheduler does not know the size of the task in advance, it is natural to take decisions as a function of w (the work already executed) instead of as a function of time, as done in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF]. To see this, let us consider Figure 1, where speed s 1 is used for the first w 1 cycles, s 2 from w 1 to w 2 cycles and so forth. The first curve shows the cumulative amount of work w(t) executed by this schedule up to time t while the second curve shows τ (w), the time needed to execute w amount of work under the same schedule. The slope of the first curve is the current speed while the slope of the second one is the inverse of the current speed (as a function of w).

t • • • •• D w • • • • • W max • • • • • w(t) w • • • • • w 1 w 2 w 3 W max t • • • • • D • • • • • τ (w)

D. Objective Function

Let us define the following quantities:

• w := w(t), the amount of work executed at time t;

• s := s(w) ∈ {s 1 , . . . , s N }, the speed used once exactly w work units have been executed, w ∈ [0, W max ]; • Q(s), the power consumption per unit of work when operating at speed s.

We have

Q(s) := P (s) dt dw = P (s) s , (2)
where derivatives are always intended as right derivatives.

Since P (0) = 0 and P is convex, we notice that Q is nondecreasing.

Our objective is to find a speed profile s : [0, W max ] → {s 1 , . . . , s N } that minimizes the mean energy consumption to execute the given job constrained by the deadline D. The mean energy consumption under speed profile s is defined by the following formula. By conditioning on the size (W = w) of the task, we get

E(s) := Wmax 0 w 0 Q(s(x)) dx dP W (w) (3)
= Wmax 0 Q(s(x)) Wmax x dP W (w) dx (4)
= Wmax 0 F c (x) Q(s(x)) dx (5)
where the second equality follows by changing the order of integration. Therefore, our objective is to solve the following infinite-dimensional optimization problem:

min s:[0,Wmax]→{s1,...,s N } Wmax 0 Q(s(w)) F c (w) dw (6)
subject to:

Wmax 0 dw s(w) ≤ D. (7)
Note that the constraint in [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF] simply states that the job has to be completed before the deadline D. To see this, let τ (w) denote the time to execute the first w epochs, as in Figure 1.

Then, dτ (w) dw = 1 s(w) [START_REF] Gaujal | A Pseudo-Linear Time Algorithm for the Optimal Discrete Speed Minimizing Energy Consumption[END_REF] and integrating both sides and using that τ (0) = 0 and that τ (W max ) ≤ D, we obtain [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF].

The constraint [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF] can be replaced by an equality without any loss of generality because a schedule that finishes before time t can always be delayed and start at time t 0 > 0 to reach D exactly; see Figure 2. This shows that such a schedule is suboptimal, as an optimal schedule will always start with active processors from time zero.

t • • • •• • D w • • • • • W max • • • • • w(t) t • • • •• D t 0 w • • • • • W max • • • • • w (t)
The problem ( 6)-( 7) is feasible if

D ≥ W max s N . (9)
Indeed, when the scheduler uses its maximal speed from the start, it must be able to complete the largest possible job before the deadline D. This condition is assumed to be true in the rest of the paper. As a side remark, notice that if the processor has no slack ( Ds N = W max ), then the solution is trivial: use all processors from the start, as long as the probability P(W ≥ W max -) > 0 , for all > 0. Dually, if D is so large that using a single processor is enough to always satisfy the constraint (D > W max /s 1 ), then the solution is also trivial: use a single processor from the start. Thus, we also assume that D ≤ W max /s 1 to make sure the constraint ( 7) is satisfied as an equality under the optimal solution.

III. A CONVEX PROGRAMMING SOLUTION

Since the scheduler has a finite set of possibilities at each w, its decisions split [0, W max ] into a set of intervals such that in each interval, the speed is constant. In general, the same number n of processors can be used in several intervals, but we will show in Proposition 2 that the optimal schedule uses an increasing number of processors with w. Therefore, an optimal schedule can be seen as a set of N -1 decisions points

x 0 = 0 ≤ x 1 ≤ x 2 ≤ • • • ≤ x N -1 ≤ x N = W max such that n processors are used in the interval [x n-1 , x n ].
With this property, the optimization of interest in ( 6)-( 7) can be transformed into the following optimization problem:

min xi, i=1,...,N -1 N i=1 Q(s i ) xi xi-1 F c (w)dw (10a)
subject to:

N i=1 x i -x i-1 s i ≤ D (10b) x i ≤ x i+1 , ∀i = 0, . . . , N -1. ( 10c
)
where x 0 = 0 and x N = W max . This optimization yields speed profiles where at most N -1 changes of speed are allowed. Specifically, the i-th change is performed when the amount of executed work reaches x i , and within the interval [x i-1 , x i ), the speed is s i . Thus, if x i-1 = x i , then speed s i is never used.

Proposition 1. The optimization problem (10) is convex, with a differentiable objective function when F is continuous. In addition, if F is strictly increasing, then it is strictly convex.

Proof. Let F := F c (w)dw. Note that the objective in (10a) can be written as

Q(s N ) F (x N ) -Q(s 1 ) F (x 0 ) + N -1 i=1 (Q(s i ) -Q(s i+1 )) F (x i ) (11)
Here, we note that F (w) is concave because F c (w) is positive non-increasing. Within the assumptions on the s n 's, Q(s i ) < Q(s i+1 ), and given that the sum of convex functions is convex, the optimization (10a)-(10c) is convex. In addition, if F is strictly increasing, then F (w) is strictly concave and the optimization (10a)-(10c) is strictly convex. Finally, the differentiability of the objective function follows by the fact that F is continuous.

Therefore, an optimizer of ( 10) can be computed efficiently by using existing algorithms from optimization theory such as the interior point method, which are guaranteed to converge to an optimal solution; see, e.g., [START_REF] Boyd | Convex Optimization[END_REF].

The next theorem, our main result, connects the optimization [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF] with the original optimization ( 6)- [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF]. It shows that any speed profile induced by the former also solves the latter.

Theorem 1. Let s * be a speed profile such that

s * (w) = s i if and only if w ∈ [x * i-1 , x * i ) (12)
where x * solves [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF], for all i = 1, . . . , N . Then, E(s * ) is an optimal schedule.

Proof. First, we show in the following proposition that we can restrict to non-decreasing speed profiles, i.e., s * (w) is nondecreasing.

Proposition 2. Any optimizer s of (6)-( 7) is non-decreasing.

Proof. We show that we can always improve a profile that is not non-decreasing. This means that, if we denote by s a feasible speed profile which is not non-decreasing, we can always find another feasible speed profile s different from s such that E(s ) ≤ E(s).

Let s be a feasible speed profile such that s(w) = s k when w ∈ [w a , w b ) and s(w) = s h when w ∈ [w b , w c ), for some w a , w b and w c such that 0 ≤ w a < w b < w c ≤ W max and for some k ∈ {2, . . . , N } and h < k. Now, consider a new feasible speed profile s that is identical to s outside the interval [w a , w c ) and such that s (w) = s h when w ∈ [w a , w b ) and s (w) = s k when w ∈ [w b , w c ), for some w b . Note that there exists a unique choice for w b and necessarily

w b < w b ; see Figure III. • w a • w b • w b • w c w t • 1 s h • 1 s k • τ (w) • 1 s k • 1 s h • τ (w)
be the random variable denoting the energy consumption under speed profile s for a task of size W . We have the following cases

• If W < w a , then E(s) = E(s ); • If W ∈ [w a , w b ), then Q(s k ) > Q(s h ) implies E(s) = wa 0 Q(s(w)) dw + Q(s k )(W -w a ) > wa 0 Q(s(w)) dw + Q(s h )(W -w a ) = E(s ) • If W ∈ [w b , w b ), then Q(s k ) > Q(s h ) implies E(s) = wa 0 Q(s(w)) dw + Q(s k )(W -w a ) > wa 0 Q(s(w)) dw + Q(s h )(w b -w a ) + Q(s k )(W -w b ) = E(s ). • If W ∈ [w b , w c ), then E(s) = wa 0 Q(s(w)) dw + Q(s k )(w b -w a ) + Q(s k )(w b -w b ) + Q(s h )(W -w b ) > wa 0 Q(s(w)) dw + Q(s h )(w b -w a ) + Q(s k )(w b -w b ) + Q(s k )(W -w b ) = E(s )
where the inequality follows because w b -w a = w c -

w b > W -w b . • If W > w c
, then E(s) = E(s ) because the time to execute the first w c epochs is identical within s and s , and because s(w) = s (w) when w ≥ w c . Since E(s) ≥ E(s ), then s cannot optimize ( 6)- [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF].

The main consequence of Proposition 2 is that we can restrict to speed profiles that change speed at most N -1 times. Thus, let x i ∈ [0, W max ] for all i = 1, . . . , N -1 be the points of the task size where the speed is allowed to change, with x i ≤ x i+1 . Let also x 0 := 0 and x N := W max . The interpretation is that the work in the interval [x i-1 , x i ) is performed with speed s i , for all i = 1, . . . , N . Within this setting, s(w) is a cadlag piece-wise constant function and

E(s) = N i=1 Q(s i ) xi xi-1 F c (w) dw (14)
and the constraint (7) can be written as

N i=1 x i -x i-1 s i = D. (15)
Therefore, the optimization in ( 6)-( 7) boils down to the optimization problem [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF]. As mentioned before, in (10b), equality has been replaced by "≤" as this does not change the set of optimizers.

IV. EFFICIENT COMPUTATION

Theorem 1 shows that solving ( 6)-( 7) amounts to solving the finite dimensional (N -1 variables) convex optimization problem in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF], and this can be done using existing algorithms from optimization theory. We notice that if the task size distribution F is continuous, then the objective function is differentiable and an optimizer can be computed by using, e.g., interior point methods [START_REF] Boyd | Convex Optimization[END_REF]. However, if the task size distribution F is discontinuous, then the objective function is not differentiable. In this case, one can always use a linear or nonlinear interpolation to obtain a continuous version of F that approximates the given F arbitrarily well. Otherwise, one can rely on existing algorithms from derivative-free optimization theory [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF].

A. Structural Properties and Algorithm

We now investigate structural properties possessed by an optimal speed profile when the objective is strictly convex and regular. These will lead to the design of a low complexity algorithm for the computation of the optimal solution.

The first result says that only consecutive speeds are used. Specifically, if the job starts under speed, e.g., s 2 and then moves to speed s 3 , either it runs with speed s 3 until it terminates or it moves to s 4 , but there is no possibility to move directly from s 3 to any s i with i ≥ 5.

Proposition 3. Assume that F is continuous and P is strictly convex. Let s * be an optimal schedule. Then, s * uses a consecutive set of speeds.

Proof. We assume that x * i-2 < x * i-1 = x * i < x i+1 for some i ∈ {2, . . . , N }, i.e., speed s i is never used while speeds s i-1 and s i+1 are both used, and show that we can improve the objective function. Let x coincide with x * except on coordinates i and i -1 where we let x i-1 = x * i -ε and

x i = x * i + q(ε). Here, ε > 0 is such that x * i -ε > x * i-2 , x * i + q(ε) < x * i+1
and another condition to be set later (see ( 16)), and q(ε) ensures that times needed to complete x * i + q(ε) units of work under the speed profiles induced by x and x * are identical. A simple computation yields (see Figure 4).

q(ε) = ε s i s i-1 s i+1 -s i-1 s i+1 -s i -1 . Outside the interval [x * i -ε, x * i +q(ε)], the objective function (10a) is identical in both x and x * . When restricted to [x * i - ε, x * i + q(ε)], the objective (10a) under x resp. x * is O x := Q(s i )
x * i +q(ε)

x * i -ε F c (w)dw O x * := Q(s i-1 ) x * i x * i -ε F c (w)dw + Q(s i+1 )
x * i +q(ε) By continuity of F c ,

x * i F c (w)dw. • • ε q(ε) • w t • 1 si-1 • 1 si+1 • τ (x * ) 1 si τ (x)
x * i +q(ε) x * i -ε F c (w)dw = (ε + q(ε))F c (x * i ) + r 0 (ε) x * i x * i -ε F c (w)dw = εF c (x * i ) + r 1 (ε)
x * i +q(ε)

x * i F c (w)dw = q(ε)F c (x * i ) + r 2 (ε),
where r 1 (ε), r 2 (ε) and r 0 (ε) are all negligible w.r.t. ε, meaning r a (ε)/ε goes to 0 when ε goes to 0, for a ∈ {0, 1, 2}.

Using the value of q(ε) and setting r(ε)

:= r 0 (ε) -r 1 (ε) - r 2 (ε), 1 F c (x * i ) (O x -O x * -r(ε)) = ε(Q(s i ) -Q(s i-1 )) -q(ε)(Q(s i+1 ) -Q(s i )), = ε P (s i ) s i-1 1 + s i -s i-1 s i+1 -s i - P (s i-1 ) s i-1 - P (s i+1 ) s i-1 s i -s i-1 s i+1 -s i = ε s i-1 (s i+1 -s i ) (s i+1 -s i-1 )P (s i ) -(s i+1 -s i ) P (s i-1 ) -(s i -s i-1 ) P (s i+1 ) . Since P is strictly convex, (s i+1 -s i-1 ) P (s i ) - (s i+1 -s i ) P (s i-1 ) -(s i -s i-1 P (s i+1 )) < 0. Therefore, O x -O x * = -εγ + r(ε)
, where γ > 0 does not depend on ε.

As mentioned before, r(ε)/ε goes to 0 when ε goes to 0. Therefore, as soon as ε is small enough so that

r 0 (ε) -r 1 (ε) -r 2 (ε) < εγ, (16)
then O x -O x * < 0.

Proposition 4. Assume that F is continuous and P is strictly convex. Then, the optimal schedule always uses the maximal speed v N .

Proof. By contradiction, assume that the last speed used in the optimal schedule is s N -1 , over the interval [x N -1 , W max ]. Now, let us consider an alternative schedule that only changes the speed in [x N -1 , W max ] and uses speed s N -2 over

[x N -1 , z 1 ], speed s N -1 over [z 1 , z 2 ]
and speed s N -2 over [z 2 , W max ] and achieves the same total work. This is always possible if z 1 is close enough to x N -1 , since s N -2 < s N -1 < s N . Using the fact that (by the deadline constraint)

z 1 -x N -1 s 1 + z 1 -x N -1 s 1 + z 1 -x N -1 s 1 = W max -x N -1 s 2 , z 2 = W max -x N -1 -s N (s N -1 -s N -2 ) s N -2 (s N -s N -1 ) z 1 . Now, the average cost over the interval [x N -1 , W max ] of the new schedule is Q(s N -1 ) z1 x N -1 F c (w)dw + Q(s N -1 ) z1 x N -1 F c (w)dw + Q(s N -1 ) z1 x N -1 F c (w)dw.
Since F c is continuous, then the average cost is differentiable. Its derivative with respect to z 1 is

F c (z 1 )(Q(v N -2 ) -Q(v N -1 ) + F c (z 2 )(Q(s N ) -Q(v N -1 ) s N (s N -1 -s N -2 ) s N -2 (s N -s N -1
) .

When z 1 = 0, the derivative is non-positive because F c (z 2 ) = 0 and strictly negative when P is strictly convex, because this implies that Q is strictly increasing. In turn, this implies that the average cost is not minimal at z 1 = 0 and that the new schedule is better than the original one for some z 1 > 0, also implying that speed s N is used by the optimal schedule. Now, given m ∈ {1, . . . , N }, let us construct the speed profile v := v m : [0, W max ] → {s m , . . . , s N } as follows:

v(w) = s i if and only if w ∈ [y i-1 , y i ) (17)
for all i = m, . . . , N where y 0 = • • • = y m-1 = 0, y N = W max and the speed change vector point (y m , . . . , y N -1 ) ∈ [0, W max ] N -m satisfies

y i = (F c ) -1 Γ m Γ i F c (y m ) , i = m + 1, . . . , N -1 (18)
where (F c ) -1 (y) := sup{z : F c (z) = y},

Γ i := s i P (s i+1 ) -s i+1 P (s i ) s i+1 -s i (19)
and

N i=m y i -y i-1 s i = D. (20)
The next proposition gives a property on v.

Proposition 5. Assume that F is continuous and P is strictly convex. Then, there exists a unique optimizer that solves (6)- [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF]. In addition, if such optimizer uses all the speeds s 1 , . . . , s N , then it is given by v with m = 1.

Proof. Uniqueness is trivial because ( 6)-( 7) is a strictly convex optimization problem. The proof is a direct application of the KKT conditions to our optimization problem. By introducing the multipliers ν ∈ R M + and λ ∈ R + , for all i = 1, . . . , N -1, the Lagrangian function associated to ( 6)-( 7) is

L(x, λ, ν) := N i=1 Q(s i ) xi xi-1 F c (w)dw + λ N i=1 x i -x i-1 s i -D + N -1 i=0 ν i (x i -x i+1 ).
Now, we obtain ∂L ∂xi = 0 if and only if

(Q(s i ) -Q(s i+1 ))F c (x i ) + λ 1 s i - 1 s i+1 + ν i -ν i-1 = 0.
Assuming that all speeds are used in the optimal solution, necessarily ν i = 0 for all i by complementary slackness and the unique optimizer must satisfy

(Q(s i+1 ) -Q(s i ))F c (x i ) = λ 1 s i - 1 s i+1 . ( 21
)
We obtain that the optimizer of (10a)-( 10b)-( 10c) is given by ( 18) and (20), as desired.

By combining the previous propositions, we can design an algorithm that computes the optimal schedule even when some speeds are not used in the optimal solution; see Algorithm 1. Specifically, we have the following result.

Algorithm 1: Dichotomy over the set of speeds for the computation of the optimal schedule.

Input: The set of speeds s 1 , . . . , s N , the power function P , the deadline D and the probability distribution of the task size F . Output: The optimal subset of speeds s U , . . . , s N , and schedule y U , . . . y N . 1 U := N ; L := 1; 2 while U > L do Theorem 2. Assume that F is continuous and P is strictly convex. Then, Algorithm 1 computes the (unique) optimal schedule. The complexity of this algorithm is log 2 (N ) times the complexity of solving the one-dimensional equation (20).

Proof. One of the technical difficulties with (18) is that it is not defined when y 1 < 0 or y 1 > W max . To deal with this, let us extend F c by any continuous strictly decreasing function on R and let us consider the solution of ( 18)-(20) using speeds s m , . . . , s N . Now, let us notice that since Γ i are strictly increasing in i and (F c ) -1 (x) ≥ {z : F c (z) ≤ x}, then if 0 < y m < W max , then 0 < y m < y m+1 < • • • < y N < W max . Therefore, if 0 < y m < W max , then the optimal schedule with speeds s m , . . . , s N is the solution y m , . . . , y N . In this case, adding lower speeds can lead to a better schedule in terms of energy. In Algorithm 1, this case is treated in Lines 7-8: the lowest speed is lowered from y m to y m with m := (L+m)/2 . If y m > W max then W max /s m < D. This implies that using the current smallest speed y m is enough to complete the task before its deadline and this is obviously the optimal solution with the current set of speeds. In this case as well, adding lower speeds can also lead to a better schedule in terms of energy. In Algorithm 1, this case is also treated in Lines 7-8.

Finally, if y m ≤ 0, there does not exist an optimal solution where the speeds s m . . . s N are all used. Since speed s N is always used in the optimal solution (Proposition 4), one must remove the low speeds to get a feasible schedule and change m to m = (U + m)/2 (Lines 5-6 in Algorithm 1).

Let us comment on the computational complexity of Algorithm 1. First, it is clear that it halts after no more than log 2 N iterations. Then, for each iteration, a one-dimensional equation needs to be solved, i.e., (20). Rearranging terms and using (18), this equation can be rewritten as

D = W max s N + N -1 i=m y i S i = W max s N + N -1 i=m (F c ) -1 Γ m Γ i F c (y m ) S i (22)
where S i := 1 si -1 si+1 . We notice that (22) admits a unique solution and that F is differentiable almost everywhere because it is increasing. If derivatives can be computed, one can solve (22) by using (efficient) standard root finding algorithms such as the Newton-Raphson or the secant methods.

B. Explicit Calculation of Speed Change Points

Within specific choices of the task size distribution F , an expression for the solution of ( 18)-(20), and thus of an optimal speed profile, can be constructed explicitly. For illustration purposes, let us consider the case where F takes the form

F (w) = 1 -K 1 - w W max q , q ∈ R where K := 1 -Wmin Wmax -q is a normalizing constant. If q = 1, we notice that W is uniformly distributed over [W min , W max ]. Then, (F c ) -1 (y) = W max -W max y K 1/q
and substituting in (18), we obtain, for i = m + 1, . . . , N -1

y i = W max -W max 1 K Γ m Γ i F c (y m ) 1/q = W max -W max Γ m Γ i 1/q 1 - y m W max .
Thus, the dependence among the y i 's is linear, and an explicit formula for y m can be obtained by using (20), which in this case boils down to a linear equation with one unknown. Using ( 22)

D = W max N -1 i=m 1 - Γ m Γ i 1/q 1 - y m W max S i + W max s N which gives y m = W max -W max 1 s N -D Wmax + N -1 i=m S i N -1 i=m S i Γ1 Γi 1/q . V. GENERALIZATIONS
We discuss some possible generalizations of our model and results, though a detailed analysis is left as future work:

• Delayed Processor Activation. It may be the case that the activation of a processor requires some setup time. Thus, if the decision of activating a processor is taken at time t, then that processor can actually be used from time t + δ, where we may suppose that the setup time δ is deterministic. In this case, a possible solution would consist in considering the speed profile s * proposed in Theorem 1 and activating a processor at time t * i -δ, where t * i is the time corresponding to the task size decision point x * i , i.e., t * i = i j=1 1 sj (x * j -x * j-1 ). This is not optimal because the processor consumes energy in [t * i -δ, t * i ) but we claim that this heuristic should work reasonably well.

• Containers. Within applications that run on top of Kubernetes, processors (or servers) run in "containers". In this setting, the activation of a single processor may trigger the activation of a number of other processors, which may remain idle but ready to use. In Proposition 3, we have shown that the optimal speed profile s * uses a consecutive set of speeds. Thus, upon activation of a new container, only one extra processor should be used if the idle processors in the container do not consume energy. However, this latter assumption may not hold in practice [START_REF] Mahmoudi | Performance modeling of serverless computing platforms[END_REF], and in this case it may be convenient to use additional processors.

VI. CONCLUDING REMARKS

We have examined the execution of a single task of unknown size on top of a multi-processor system and investigated the problem of dynamically activating the number of processors that minimizes the overall electrical energy consumed by executing the task subject to a hard deadline constraint. We have shown that this problem can be formulated as a convex optimization problem and we have identified key structural properties that allowed us to design an extremely efficient algorithm for the computation of the optimal dynamic allocation of processors.

Our work is in close relationship with [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF]; see Section I-A. In this reference, we stress that the speeds available to the scheduler form an interval of the real numbers, while we consider a finite number of speeds, i.e., N . Under appropriate conditions, it would be interesting to investigate whether the proposed optimal speed profile converges to the optimal speed profile proposed in [START_REF] Lorch | Improving dynamic voltage scaling algorithms with PACE[END_REF] in the limit where N → ∞, and if convergence occurs, then at what rate. Also, our results facilitate the design of heuristics to handle the case where a stream of tasks of unknown duration join the system of over time. This case is much more difficult to deal with: the hard deadline constraint makes non-trivial even just the problem of finding a feasible schedule. We leave these natural continuations of our work for future research.

Figure 1 .

 1 Figure 1. The function w(t) is the cumulative work executed by time t under a given schedule (left) and τ (w) is the time needed to execute work w under the same schedule (right). The speed is the slope in the left figure and the inverse of the slope in the right one.

Figure 2 .

 2 Figure 2. Shifting the starting time to reach D.

Figure 3 .

 3 Figure 3. Execution time functions τ (w) and τ (w) corresponding to schedules s and s respectively. Let

Figure 4 .

 4 Figure 4. Construction of x and q(ε), starting from x * and ε.

3 m 5 if y m ≤ 0 then 6 L

 356 := (U + L)/2 ; 4 Solve (20) for y m using speeds s m , . . . , s N , where y m+1 , . . . , y N are given by (18);