Mohamed Essadki

Jonathan Jung

Adam Larat

Milan Pelletier

Vincent Perrier

A task

come

Introduction

Computations of turbulent flows has been undergoing two breakthrough over the last three decades, concerning computational architectures and numerical methods. On the one hand, computing clusters are reaching a size of many hundreds of thousands of cores, which nowadays allows to overcome RANS modeling by considering unstationary LES simulations [START_REF] Chalot | A consistent finite element approach to large eddy simulation[END_REF][START_REF] Hill | Hybrid tuned center-difference-weno method for large eddy simulations in the presence of strong shocks[END_REF]. On the other hand, this new framework requires high order numerical methods, which have also been much improved recently [START_REF] Chalot | A consistent finite element approach to large eddy simulation[END_REF][START_REF] Dumbser | A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes[END_REF]. Nevertheless, these newly manufactured massively multicore computers rely on heterogeneous architectures. In order to improve the efficiency of our methods on these architectures, their implementation and algorithmic must be redesigned. The aim of this project is to optimize the implementation of all-in-one numerical methods and models on emerging heterogeneous architectures by using runtimes schedulers. In this publication, we focus on StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF].

Runtime schedulers allow to handle heterogeneous architectures in a transparent way. For this, the algorithm must be recast in a graph of tasks. Computing kernels must then be written in different languages, in order to be executed on the different possible devices. Then, the runtime scheduler will be in charge of dynamically balancing the tasks on the different available computing resources (cores of the host, accelerators, etc.).

Finite differences methods have already been successfully ported on hybrid architectures [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to nvidia graphics cards using cuda[END_REF]. However, we aim at developping numerics on unstructured meshes, for example of the discontinuous Galerkin type, and our methods have a very different memory access pattern, which is one of the bottlenecks of the task-driven programming. In order to minimize the complexity of our presentation, gathering both the numerical procedure and its implementation within the non-classical StarPU framework, we have decided to restrict its scope to first order in space and time.

A first order finite volume solver for two dimensional conservation laws

Hyperbolic conservation laws

Let Ω be a subset of R 2 . We consider a two dimensional system of conservation laws with source terms

∂ t W + ∇ • F(W ) = S(W ), ( 1
)
where t is the time, ∇ = (∂ x , ∂ y ) T with (x, y) the spatial position, W ∈ R m (m ∈ N) is the vector of conservative variables, F = (F x , F y ) T : Ω → R m is the flux function and S(W ) represents the source terms. Later, the number of conserved variables is noted nVar = m. The unknowns W depend on the spatial position (x, y) and on the time t. We assume that system (1) is hyperbolic, meaning that the directional Jacobian ∂F ∂W • n is diagonalizable in every direction n ∈ S 1 , and we note λ 1 (W ) ≤ • • • ≤ λ m (W ) its eigenvalues, the direction n being implicit. We aim at computing a solution of (1) with initial condition ∀(x, y) ∈ Ω, W (x, y, 0) = W 0 (x, y).

(2)

Finite volume discretization

In this section, we present a simple numerical scheme which solves system (1) in space and time: (x, y, t) ∈ Ω × R + . If Nx and Ny are two given positive integers, we define the sequences x i+ 1 2 0≤i≤Nx and y j+ 1 2 0≤j≤Ny by x i+ 1 2 = a + i∆x, i = 0, • • • , Nx and y j+ 1 2 = c + j∆y, j = 0, • • • , Ny, so that Ω is discretized into Nx × Ny cells C i,j = x i- 1 2 , x i+ 1 2 × y j-1 2 , y j+ 1 2 . We also consider a sequence of times (t n ) n∈N such that t 0 = 0, which defines the time steps (∆t) n = t n+1 -t n > 0.

Hyperbolic part

Next, for any cell C i,j of the mesh, at any time t n , we look for an approximation of the average value of the solution W (x, y, t) on the cell:

W n i,j = 1 |C i,j | Ci,j
W (x, y, t n )dxdy.

(3)

By integrating system (1) over a cell C i,j , one obtains:

d dt Ci,j W (x, y, t)dxdy + Γ⊂∂Ci,j Γ F(W ) • n dS = 0,
where n = (n x , n y ) is the outward unit normal to the boundary ∂C i,j . Then, the first order explicit finite volume scheme writes [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]:

W n+1, * i,j -W n i,j = - (∆t) n ∆x F W n i,j , W n i+1,j , (1, 0) -F W n i-1,j , W n i,j , (1, 0) - (∆t) n ∆y F W n i,j , W n i,j+1 , (0, 1) -F W n i,j-1 , W n i,j , (0, 1) , (4)
where F (W k,l , W m,n , n) is a chosen numerical flux. In all our computations, the first order Lax-Friedrichs flux will be used:

F (W L , W R , n) = F(W L ) • n + F(W R ) • n 2 - σ 2 (W R -W L ) with σ = max p∈ 1,m (|λ p (W L )| , |λ p (W R )|) .
The numerical scheme ( 4) is stable provided the following CFL condition is ensured

(∆t) n × max (i,j)∈ 1,Nx × 1,Ny max p∈ 1,m |λ p (W n i,j )| ≤ min(∆x, ∆y). (5)
In practice, a constant time step (∆t) n = ∆t is set at start, and therefore we just need to check at the beginning of each iteration that ∆t satisfies inequality (5).

Source term

Since we consider only a first order numerical scheme in space and time, the source term treatment can be simply added by use of a first order operator splitting technique [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF].

Therefore, the additional source terms can possibly be taken into account by local cell-wise integration of the following ODE, which is simply discretized by mean of a forward Euler time scheme:

∂ t W = S(W ) -→ W n+1 i,j = W n+1, * i,j + ∆tS W n+1, * i,j , (6)
where

W n+1, * i,j
is the partial update given by the Finite Volume integration (4). Let us note that though the update ( 6) is simple and cheap, the computation of S W n+1, * i,j may be rather complex and computationnally expensive.

StarPU: a task scheduler

Task scheduling offers a new approach for the implementation of numerical codes, which is particularly suitable when looking for an execution on heterogeneous architectures.

The structure of the codes is divided into two main layers that we will call here bones and flesh. At first level, the code is parsed to build a tasks dependency graph that can be viewed as its skeleton. Next, one enters within this graph of tasks and starts executing the actions associated to each of these tasks. Then, the code acts on the memory layout thanks to what can be viewed as its muscles. During the execution part, one may have different choices in the tasks execution on the available computational power and this is where the scheduler plays its role: it aims at distributing the actions so as to minimize the global computational time.

StarPU tasks

In StarPU, a task is made of the following: a codelet, associated kernels and data handles.

• The codelet is the task descriptor. It contains numerous information about the task, including the number of data handles, the list of available kernels implemented to execute this task (for CPU, GPU or possibly something else. . . ) and the memory access mode for each data handle: "Read" (R), "Write" (W) or "Read and Write" (RW). • The kernels are the functions that will be executed on a dedicated architecture: CPU, GPU, etc.

A task may have the choice between different kernels implementation and it is the role of the StarPU scheduler to distribute the tasks on the available heterogeneous architectures following a certain criterion (in general minimizing the global execution time). • The data handles are the memory managers. Each data handle can be viewed as the encapsulation of a memory allocation, which allows to keep trace of the action (read or/and write) of the task kernel on the memory layout. In particular, this allows to build the task dependency graph.

Construction of the tasks dependency graph

The domain of size Nx × Ny is partitioned into NPartX × NPartY balanced parts of size NxLoc × NyLoc, see Figure 1. One task will deal with one partition at a time. In order to treat the interaction between the domains in an asynchronous manner, each domain is supplemented with copies of its border data, see Figure 2.

• Each memory allocation is associated with a data handle: one for each subdomain, four for its borders copy, called overlaps in the following, • The codelet of each task points toward a certain number of these data handles, tagged as "R" or "W" or "RW", for Read and/or Write, following the access mode exepected by the associated kernel, • Finally, the code is parsed and the tasks are submitted to StarPU. When a task has to write on a certain data handle memory, a dependency edge is drawn between this task and the latest tasks having read access on this data handle. Similarly, for each of its read-accessed data handles, a task will depend on the latest tasks having write access on them. • Once the dependency graph is built, the tasks are distributed in dependency order by the scheduler on the available computational ressources, following the available kernel implementations. Note that building the dependency graph and lauching the task can be an intricated work, especially since the future of the graph may depend on the computation itself: think of the number of time steps which may depend on the solution, for example.

Task overhead

StarPU tasks management is not free in terms of CPU time. In fact, each task execution presents a latency called "overhead ". We want to know when this additional time can be neglected or not. There is a StarPU benchmark that allows to measure the minimal duration of a task to ensure a good scalibility. We send 1000 short tasks with the same (short) duration varying between 4 and 4096 µs and we study the scalability. In Figure 3, we plot the results obtained on a 2 dodeca-core Haswell Intel Xeon E5-2680 and on a Xeon Phi KNL with the scheduler eager. On few cores, we have a good scalibility result, even if the duration of the tasks is short (< 0.2ms). On many cores, we need longer tasks duration to get satisfying scaling (≈ 1ms).

To sum up, if the duration of the tasks is smaller than the microsecond, their overhead cannot be neglected anymore.

Schedulers

The purpose of the scheduler is to launch the tasks when they become ready to be executed. In StarPU, many different scheduling policies are available. In this paper, we consider only the two following:

• the eager scheduler: it is the scheduler by default. It uses a central task queue. As soon as a worker has finished a task, the next task in the queue is submitted to this worker. • the dmda scheduler: this scheduler takes into account the performance models and the data transfer time for the tasks (see Section 5.3.2). It schedules the tasks so as to minimize their completion time by carefully choosing the optimal executing architecture.

Practical Implementation

After a general presentation of the model and the numerical method in section 2 and of the runtime environment with a particular focus on StarPU in section 3, we now detail the way we have implemented this numerical resolution of a system of hyperbolic conservation laws within the task-driven framework StarPU.

Memory allocation

Since the tasks dependency graph is mainly based on the memory dependency between successive tasks, memory allocation is a crucial development part of our application. This starts with the definition of the structure Cell, which contains the nVar = m conserved variables at a certain point in space and time.

Principal variables

For a first order finite volume resolution of a system of conservation laws, only two main variables are needed:

• u, the computed solution. In each mesh cell, it contains nDoF × nVar floats. However, since at first order nDoF = 1, u contains only one Cell structure of size nVar corresponding to one local approximation (3) of the solution per mesh cell. • RHS, the vector of residuals. It has exactly the same memory characteristics as u, since at each time step, the update is done by:

u += RHS. (7)
These two vectors of variables are allocated per subdomain. Typically, for each of the NPartX × NPartY subdomains, a vector uLoc of (NxLoc × NyLoc) Cells is created and encapsulated in a StarPU data handler which allows to follow the memory dependency.

Overlap additional memory buffers

In order to minimize the communications and the dependencies between subdomains, each subdomain comes with four additional one-dimensional buffers corresponding to the possible four overlap-data needed at the subdomain boundaries (East, North, West and South).

Two vectors of size NxLoc × sizeof(Cell) (ovlpS and ovlpN) and two vectors of size NyLoc × sizeof(Cell) (ovlpE and ovlpW) are always allocated, whether the overlap needs to be used or not. The reason for that is that the number of data handlers passed to a StarPU kernel needs to be constant. Therefore, when copying the overlap-data for example, all the overlap data handlers are passed anyway but nothing is done if they are not needed, like in the case of a periodic subdomain in one direction.

Of course, here lies a small communication optimization, when sending a task on another device, since some useless data is transfered. However, we think that the overlap tasks should be of negligeable size compared to the task acting on the entire subdomains and should be mainly executed on the host node.

Description of the tasks

In paragraph 2.2.1, we have seen that the first order finite volume method essentially consists in 1) checking the computational time step ∆t, 2) computing all the numerical fluxes at all the edges of the mesh, 3) gathering them in the RHS vector and finally 4) updating the numerical solution u.

Based on this general decomposition, here is the list of tasks implemented in our application. Between brackets is specified the memory data handlers accessed by the task, with their respective access rights given between parentheses:

• initialCondition uLoc(W) : fills each subdomain solution with the initial condition.

• checkTimeStep uLoc(R) : computes the largest characteristic speed within the subdomain. In order to avoid gathering these time constraints globally, we only check that the fixed time step ∆t initially given by the user respects locally the stability constraint. Copy to overlaps tasks, example with two domains. The global computational domain is vertically divided into two parts, blue (left) and green (right). The green one is supplemented with a west overlap, whereas the blue one is supplemented with an east overlap. One residual computation includes two communications tasks: copying the right column of the blue domain into the west overlap of the green domain, and copying the left column of the green domain into the east overlap of the blue domain.

East West

• copyOverlaps ovlpE(W),ovlpN(W),ovlpW(W),ovlpS(W),uLoc(R) : each subdomain fills the corresponding neighbors overlap data vectors, if needed. This is not the case if the subdomain is periodic in one or both directions, or if one of the boundary states is prescribed (Dirichlet boundary condition).

The northest (resp. southest) line of the subdomain is copied in the ovlpS (resp. ovlpN) data handle of the northern (resp. southern) neighbor. The extreme east (resp. extreme west) column is copied in the ovlpW (resp. ovlpE) data handle of the eastern (resp. western) neighbor. Figure 4 depicts the copyOverlaps task for two domains. • internalResiduals uLoc(R),RHS(W) : computes the numerical flux at all the edges of the subdomain, except those of the boundaries where an overlap data vector has to be used:

RHS i,j = ∆t ∆x F * i+ 1 2 ,j -F * i-1 2 ,j + ∆t ∆y F * i,j+ 1 2 -F * i,j-1 2 . (8)
• borderResiduals ovlpE(R),ovlpN(R),ovlpW(R),ovlpS(R),uLoc(R),RHS(RW) : computes the remaining numerical fluxes, meaning those between the overlap vector states and the border cells of the subdomain. Therefore, the overlap data vectors passed in argument here are those belonging to the subdomain. • update uLoc(RW),RHS(R) : update the numerical solution subdomain-wise, thanks to the update relation [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to nvidia graphics cards using cuda[END_REF]. The corresponding task diagram for one time step and two sub-domains is given in Figure 5. This graph is an output we can get from StarPU to verify the correct sequence of the tasks.

Specific asynchronous treatment of the outputs

The output task consists in writing the global mesh solution u into a disk file at once. In order to avoid global synchronization, a data handle dataForOutput of size Nx × Ny × nVar is used to temporarily store the data of each subdomain. Since this buffer will be entirely written in an output file by a single task outputTask, it needs to be contiguous in memory and to be encapsulated in a single global descriptor. However, when the solution needs to be written on the disk, it is first copied in this dataForOutput buffer and this should be done in an asynchronous parallel manner, by the gatherForOutput tasks. This is only possible if the dataForOutput buffer is also described per subdomain, see Figure 6. At the end, we need one global data handle and a set of NPartX × NPartY per-subdomain data handlers, both sharing the same memory allocation. Since the dependency graph is built on the memory dependancy between the tasks, it is very dangerous to use data handlers sharing the same memory slots.

Fortunately, both descriptors can be declared in StarPU and a specific procedure allows to switch from one to another, so that they are never used concurrently. First, the global buffer is allocated and encapsulated in a global StarPU data handler. Next, it is also described as a vector of NPartX × NPartY data handlers, thanks to the command: starpu_matrix_data_register(starpu_data_handle *data_handle,int DEVICE,void *data_ptr, int LD,int NxLoc, int NyLoc, size_t sizeof(Cell));

Here, the integer LD allows to encapsulate the grid block per block, since the global vector is accessed through the formula:

for(int j=0;j<NyLoc;j++) {for(int i=0;i<NxLoc;i++) global_data[i+j*LD];} Then, if data_ptr points to the bottom-left corner of the subdomain, the corresponding block is encapsulated in the subdomain data handler data_handle.

Eventually, when the solution needs to be written on the disk, the subdomain-wise descriptor is activated and each subdomain copies its local solution uLoc to the appropriate dataForOutput sub-data-handle (gatherForOutput task). Next, the description of the dataForOutput buffer is switched to its global descriptor and the write-to-disk single task can be executed (outputToDisk task). Once the writing is finished, the description of the memory buffer is switched back to the partitioned sub-data-handles. Figure 7 illustrates how the simple task diagram for the numerics shown in Figure 5 is updated for outputs. In Figure 8, we plot the Gantt chart obtained for 60 time iterations of the finite volume scheme on two cores with an output every 20 iterations. Green color indicates working time and red color symbolizes sleeping time. This asynchronous treatment of the outputs allows to perform outputs without blocking one core for the outputs and without inducing sleeping time for the other cores.

Computing data handlers

Keep computation

Output data handlers Unpartition Partition Figure 6. Strategy for non blocking IO. The computing buffers are depicted on the left. The data handlers for output dataForOutput are depicted on the right. Each part of the domain is writing on a dedicated part of dataForOutput using the task gatherForOutput. Once this copy is done, the computing buffers can be used for further computations. On the output handlers side, once all the copies have been completed, the buffer can be switched to a global buffer thanks to a switch task, which is followed by a starpu_data_invalidate_submit command on the local descriptor of dataForOutput. Once the unpartitioning is done, the output handler can be used as a single data handler, and one tasks is in charge of writing this buffer in an output file. Once the ouput in the file is finished, another switch task is in charge of repartitioning the data, the global descriptor is invalidated, and the dataForOutput is able to receive data from the computing buffers again.

Results

In this last section, we study the performance of our task-driven implementation on test cases solving the two-dimensional Euler equations over the periodic domain R Z where E is the total energy, p is the pressure and H = E + p/ρ is the total enthalpy. In the following, we assume that p satisfies a perfect gas pressure law p = (γ -1)ρ(E -u 2 +v 2

2

). The source term is set to zero: S(W ) = 0. This system of equations ( 9) is known to be hyperbolic with characteristic velocities given by u.n -c, u.n, u.n and u.n + c, c being the sound velocity given by c 2 = γp/ρ.

Test cases

On this set of equations, we use two classical test cases, namely a linear advection of a density perturbation and an isentropic vortex. Both have the advantage to present an analytical solution.

Localized cosine perturbation on ρ

First, we consider an arbitrary perturbation of the density field ρ(x, y) which is simply advected on a constant field ū, v, p. We choose: r(x, y) = (x -0.5) 2 + (y -0.5) 2 , ρ(x, y) = ρ (r(x, y)) = 1 + (r(x, y) < 0.25) * cos (4πr(x, y)) , ū = v = 1, and p = 1/γ, [START_REF] Spiegel | A survey of the isentropic euler vortex problem using high-order methods[END_REF] so that the sound velocity is everywhere smaller than one.

At any computational time t, the exact solution is given by r(x, y, t) = (x -ūt -x -ūt -0.5) 2 + (y -ūt -y -ūt -0.5) 2 , ρ(x, y, t) = ρ(r(x, y, t)), [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes[END_REF] where the operator • stands for the floor function, in order to take the periodicity into account.

Isentropic vortex

Next, we consider a slightly more complex test case which couples all four equations but still provides a quasi-analytical solution. The "quasi " prefix is explained later on. The initial solution is isentropic and its velocity field is a superposition of a constant and a divergence-free fields:

s = p ρ γ = 1 γ , u(x, y) = ū + ũ(x, y), s.t. ∇ • ũ = 0. ( 12
)
From these two last statements, it is obvious that the entropy and the density are simply advected:

D t s = 0 and D t ρ = 0, (13)
D t . = ∂ t . + u • ∇.
being the material derivative. Now, one would like the velocity field to be simply advected by its constant part ū, and this comes if and only if:

(ũ • ∇) u + ∇p ρ = 0. (14)
This is verified on R 2 by an isentropic rotating Gaussian field:

ω(r(x, y)) = ω exp - r 2 (x, y) 2R 2 , ũ(x, y) = - y R ω(r), x R ω(r) t , ρ(r) = 1 - γ -1 2 ω 2 (r) 1 γ-1 , ( 15
)
where r is the same radius function as in [START_REF] Spiegel | A survey of the isentropic euler vortex problem using high-order methods[END_REF], R is a characteristic radius of the central perturbation and ω is the vortex intensity. The numerical values of our setup are:

ω = 1, R = 0.1, ū = v = 1, γ = 1.4. (16)
Since the density field is radial, the divergence-free field does not modify it and the whole initial solution is just advected by the constant velocity field ū = (ū, v) t , so that the analytical relation [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes[END_REF] almost still stands true. In fact, the whole previous development is true over R 2 but the connection of the velocity field at the boundaries of our periodic domain Ω is not continuous. Nonetheless, the jump exponentially decreases when the size of the domain increases, which comes exactly to the same as decreasing R. Then, for a given final time T , we can choose R such that the error of our approximated analytical solution is negligible compared to the numerical error, so that the latter can be measured, [START_REF] Shu | Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws[END_REF][START_REF] Spiegel | A survey of the isentropic euler vortex problem using high-order methods[END_REF].

Validation by mesh convergence study

In Figure 9 we show the convergence results obtained on both test cases described above. This validates the implementation of our numerical method since it is of order 1, as expected.

Parallel efficiency

We perform a strong scaling study. It means that we keep the same problem size and we increase the number of cores. In Figure 10, we consider a mesh of 1024 × 1024 cells subdivided in NPart = 1, 4, ..., 16384 sub-domains and we perform a scalability study on two different architectures: a 2 dodeca-core Haswell Intel Xeon E5-2680 and a Xeon Phi KNL. Recall that the number of tasks is proportional to the number of sub-domains NPart.

To have a correct scaling, we need enough tasks. NPart should be at least of the order of the number of cores. However, if we increase the number of tasks NPart a lot, the size (time duration) of each task becomes to small compared to the task overhead (see section 3.3) and the scaling starts to saturate. To sum up, in order scale correctly on many cores, we need a lot of tasks of long duration (compared to the tasks management overhead) and therefore, we need big meshes.

In Figure 10, we distinguish three types of curves. Those for which we do not provide enough tasks: NPart = 1, 4, 16 in 10a and 10b. Those which scale correctly on all the cores: NPart = 64, 256 in 10a, NPart = 64 in 10b. And finally, those who start to saturate due to task overhead: NPart = 1024, 4096, 16384 in 10a and the same partitions including NPart = 256 in 10b.

GPU implementation and performance

GPU implementation

In order to launch part of the computation on GPUs, one just has to write a CUDA or OpenCL version of the kernels, see paragraph 3.1. We choose to use CUDA. Then, we let the scheduler choose between a CPU device or a GPU device, for each task. It is not necessary to have a CUDA implementation of each kernel, but for performance reason (to minimize the communication between the CPU and the GPU) we need at least to write a CUDA version of the following kernels: checkTimeStep, copyOverlaps, internalResidual, borderResiduals and update, see paragraph 4.2. Let us detail the CUDA implementation of these kernels:

• copyOverlaps: we associate a virtual processor for each variable of each cell and we copy uLoc into olvpE, olvpN, olvpW or olvpS. For olvpN and olvpS, since two successive processing elements access two neighboring memories (in global memory), it allows to achieve optimal memory bandwidth (coalescing access). This is not the case for olvpE and olvpW. • update: we associate a virtual processor for each variable of each cell and we add RHS to u, see [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to nvidia graphics cards using cuda[END_REF]. The memory accesses are also coalescent. • checkTimeStep: we associate a virtual processor for each cell. The checking is a little bit different from the CPU implementation. Indeed, computing the largest characteristic speed is quite difficult in parallel. Then, we compute a local time step

(∆t) i,j = min(∆x, ∆y) max p∈ 1,m |λ p (W i,j )|
and we check if (∆t) i,j is greater than the constant time step ∆t. • internalResidual and borderResidual: we associate a virtual processor for each cell. We compute the numerical flux on each edge of this cell and we add the result to the right hand side, RHS, of this cell. In the GPU implementation, we compute the fluxes twice: the two vitual processors associated to cells C i,j and C i+1,j compute the same flux F * i+ 1 2 ,j . Indeed, computing the flux once at each edge and balancing it on the two neighboring cells would induce concurrent memory access. For example, the fluxes computed at edges Γ i- 1 2 ,j and Γ i+ 1 2 ,j could be added simultaneously to RHS i,j . However, since computation is fast on GPU, even if we compute all the fluxes twice, we still observe good performance.

Performance model

To reach good scheduling performances on heterogeneous architectures (with CPU and GPU for example), StarPU needs to be able to estimate in advance the duration of a task. This is done by associating each codelet with a performance model. This performance model provides for each codelet and for a certain range of data sizes, an expected average and standard deviation of the task execution times.

We use this tool to compare the execution time of different tasks on a CPU Haswell Intel Xeon E5-2680 and on a GPU Nvidia K40-M. We consider the numerical test of Section 5.1.1. The computational domain is divided into 8192 × 8192 cells. In Figures 1112, we plot the average execution times on CPU and GPU of the tasks associated to a specified codelet as a function of the tasks data size. To decrease the tasks data size, we increase the number of sub-domains from 2 × 2 up to 256 × 256. The maximum size of the tasks is fixed by the memory of the GPU: 10 GB in our case.

From Figures 1112, we immediately notice that computationally expensive tasks, such as checkTimeStep, internalResidual and update, are greatly accelerated on GPU. Observed speedups are between 2.6 and 24 for internalResidual and between 0.33 and 20 for update. However, tasks involving more data transfer than useful calculations, such as every task associated with subdomains boundaries (copyOverlap and borderResidual), present no reason to be executed on a GPU in the range of data size commonly used. This is where the scheduler starts playing an important role: the overall computation will accelerate if part of its tasks are executed on the GPUs, but certainly not all of them.

Graph of activities with CPUs and GPUs

With StarPU, there are different tools to study the activities of each worker. We look here at the graph of activities: it shows the activity of each worker and the evolution of the number of tasks available on the system with respect to the execution time. In Figure 13, we plot the graph of activities obtained with a run of 50 iterations of the finite volume scheme on a 16384 × 16384 cells domain, subdivided in 16 × 16 sub-domains, on a heterogeneous architecture of 4 CPUs and 4 GPUs. We consider the work done by two different schedulers: eager and dmda. Green sections indicate time spent on the kernels execution. Red sections indicate the proportion of time spent in StarPU. Black sections indicate sleeping time. With both schedulers, our graphs are essentially green, which means that all our workers have enough work: the granularity of the tasks is good and thus the parallelism is also good. We have a lot of tasks (16 × 16 sub-domains) and each sub-domain stays big (more than four millions of unknowns). However, even if we observe only little waiting time on both graphs, we remark that the execution time is more than two times smaller with the dmda scheduler than with the eager one (61 s versus 146 s). In fact, the dmda scheduler takes the performance models into account and therefore executes every task on its optimized device. These graphs are very useful to know the activities of each worker but do not give any information about the tasks executed on each worker. To obtain such information one would need to plot a Gantt chart, which is an other performance tool provided by StarPU. This is now the main content of the next section which should be considered as a concluding and opening section.

Opening: toward HPC of complex multiphase flows on large scale heterogeneous clusters

As we have seen during the last result section, the performance of the task driven implementation strongly depends on the work load of each task. In a nutshell, we need a large number of big tasks. This is simply achieved by increasing the local workflow, which could be achieved by increasing the number nVar of equations, by increasing the number nDoF of degrees of freedom per cell, or by adding some local work by mean, for example, of a source term. We start by considering the latter case.

A model for the dynamics of an evaporating disperse phase

Without entering some complex explanation, we consider the following hyperbolic system of conservation laws, [START_REF] Essadki | High order moment model for polydisperse evaporating sprays towards interfacial geometry description[END_REF]: where m k/2 (t, x) = 1 0 S k/2 n(t, x, S)dS are the fractional moments of a certain size distribution n(t, x, S), S indexing the size, evaporating at rate K and u(t, x) is its velocity field which relaxes toward an underlying gas field u g at time rate θS, thanks to the bottom right term.

           ∂ t m 0/2 + ∂ x • m 0/2 u = -K n(t, x, S = 0), ∂ t m 1/2 + ∂ x • m 1/2 u = -K 2 m -1/2 , ∂ t m 2/2 + ∂ x • m 2/2 u = -K m 0/2 , ∂ t m 3/2 + ∂ x • m 3/2 u = -3K 2 m 1/2 , ∂ t m 2/2 u + ∂ x • m 2/2 u ⊗ u = -K m 0/2 u + m 0/2 ug-u θ , (17)
✶ ✶ ✶ ✶ ✶ ✶ ✶✁✂ ✄ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ s☎✆✝✞✟✟✁✠ r✁✡✠☛ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ ✷ ✹ ✄ ✽ ✶ ✶ ✷ ✸ ✹ ✺ ✄ ✼ ✽ (
This system has now nVar = 6 equations. The left hand part is a simple linear transport along the velocity field u(t, x). The right hand side is made of two source terms: an evaportation term acting on all the moments and a drag term acting on the velocity components only, modeling the traction of the spray by an underlying gas field. As discussed in 2.2.2, the system terms are split so that the resolution of (17) boils down to

∂ t M + ∂ x .F (M ) = 0, d t M = -K [n(t, x, S)] and d t (m 2/2 u) = m 0/2 u g -u θ . (18)

Integration strategy

Following what has been said in previous sections, our numerical procedure is rather simply adapted to this new system. The only key point lies in the evaporation term which requires the reconstruction of the entire size distribution n from the moments m k/2 , k = 0, . . . , 3 by an entropy maximization procedure, see [START_REF] Essadki | High order moment model for polydisperse evaporating sprays towards interfacial geometry description[END_REF] for more details. This allows to estimate the unknown quantities n(t, x, S = 0) and m -1/2 numerical procedure being convex state preserving, the set of moments is everywhere realizable and this local reconstruction is always possible in every cell. However, the local reconstruction of n involves the computation of exponential of matrices which are rather costful, especially when compared to the simplicity of the other terms.

This explains why we have tried to pass this reconstruction procedure on GPUs and compared the computational times when activating the accelerating device or not and when using one scheduler (eager) or another (dmda). The results are presented in Table 1 and in Figure 14. All simulations hereafter compute 100 time iterations on a 200 × 200 domain, divided into 2 × 2, 2 × 4, 4 × 4 or 4 × 8 subdomains, of an evaporating spray within a Taylor-Green vortices velocity field for the gas, the same as the latest test-case of [START_REF] Essadki | High order moment model for polydisperse evaporating sprays towards interfacial geometry description[END_REF]. In table, 1, we see first that activating the GPU immediately improves the computational time a lot: the reconstruction procedure requires fast linear algebra GPUs are very good at. Next, the eager scheduler is not smart and does not anticipate the global computational time of the tasks. When switching to the dmda scheduler, we gain a little bit more performance by better distributing the tasks on the devices: this is illustrated on the Gantt charts in Figure 14. On the top figure, we look at the task scheduling when no GPU is activated. This is the eager scheduler. dmda scheduler does only worsen the situation, since each of the four CPUs is permanently loaded with tasks. We do not see it very much on this image due to the large number of tasks, but the source term tasks are about four time larger than the others: 80% of the effort is spent on these tasks. However, if we activate the GPU kernel for the source term tasks, we immediately reduce the computational time by a factor 2.6, when the eager scheduler is kept, see Figure 14b. But there, we also see that not all the source term tasks have been given to the GPU, even though it is much faster at treating these tasks. This explains why, when turning the dmda scheduler on, Figure 14c, the computational time is again reduced by more than a factor 2: now all the source term tasks are sent to the GPU and only the transport part remains on the CPU.

Nonetheless, we also see that some red stripes corresponding to slipping time, remain on the CPU line, meaning that one could gain even a little more by activating the computation of the transport tasks on the GPU card.

Conclusion and opening

We acknowledge the last results of this section are poor in term of physical meaning, especially since the numerical method is limited at first order and the model integration is not fully assessed. Nonetheless, this paper has stayed technical on purpose, in order to explore all the possibilities and restrictions of the implementation of a numerical method within the framework of a runtime, namely StarPU.

We have described and shown experimentally the main features of our code: the numerical procedure described in section 2 is subdivided into tasks (section 4), after a quick presentation of the StarPU environnement in section 3. In section 5, we have clearly demonstrated that the task distribution is efficient if the tasks are both large and numerous enough, what requires a large enough problem at start. Then, we have looked at the task distribution on a heteregeneous architecture and shown that all the tasks do not deserve to be executed on the GPU accelerator: only the one with the highest arithmetic intensity should be loaded there, while memory consuming tasks should stay on the host node. This demonstrates the necessity for an efficient scheduler which is going to make the task distribution decisions on the fly.

In this last section, we have had the will to go toward a more HPC oriented application. Instead of looking at standard Euler equations, we have implemented a model for the dynamics of an evoporating spray including intensive source terms computation. This test case is very promising for our study. Indeed, we intend to go to higher order numerical methods, which imply an increased number of degrees of freedom per cell and an increased local arithmetic intensity. Moreover, this considered model is part of a hierarchy of models which can be enriched at will by involving an increasing amount of meaningful variables and associated equations. All of this goes in the right direction, which is a larger number of computationally expensive tasks.

The work at CEMRACS16 within the Hodins team has been fruitful, with many interesting results and a common mastering of the task-driven programming for physical applications. As it is now clear, the results are promising enough for the collaboration to go on and further studies now need to be lead with higher order numerical procedures, which are going to be of the discontinuous Galerkin type.

Figure 1 .

 1 Figure 1. Partitioning of an initial mesh of {Nx = 30 × Ny = 30} cells into {NPartX = 2 × NPartY = 2} 225 cells domains.

Figure 2 .

 2 Figure 2. Partition (in blue) with its overlap (in red) in the East, North, West and South direction. Data of each partition is composed of these five handles.

  2 dodeca-core Haswell Intel Xeon E5-2680.

  Xeon Phi KNL.

Figure 3 .

 3 Figure 3. Task size overhead with the eager scheduler: scalability results obtained with duration of tasks varying between 4 and 4096 µs on two different machines.

Figure 4 .

 4 Figure 4.Copy to overlaps tasks, example with two domains. The global computational domain is vertically divided into two parts, blue (left) and green (right). The green one is supplemented with a west overlap, whereas the blue one is supplemented with an east overlap. One residual computation includes two communications tasks: copying the right column of the blue domain into the west overlap of the green domain, and copying the left column of the green domain into the east overlap of the blue domain.

Figure 5 .

 5 Figure 5. Task diagram built by StarPU for one Forward Euler step of our first order finite volume task-driven implementation on two subdomains, hence the horizontal symmetry.

2 .Figure 7 .Figure 8 .

 278 Figure 7. Task diagram built by StarPU for one time step, including output of the initial and final solutions. The solution u is first copied subdomain per subdomain into a global buffer dataForOutput, after which this buffer is written in a file by a single task. Since two different descriptors of the buffer are needed, a switch task is launched between both actions.

  (a) Localized cosine. Straight lines are estimated slopes of the order of 0.95. (b) Isentropic vortex. Straight lines are estimated slopes of the order of 0.95.

Figure 9 .

 9 Figure 9. Mesh convergence study for both test cases. Errors are L 1 (red), L 2 (blue)

Figure 10 .

 10 Figure 10. Strong scaling with the eager scheduler: scalability results obtained on two different architectures with a number of subdomain NPart varying from 4 to 16384.

  Codelet update.

Figure 11 .

 11 Figure 11. Average of the execution times on CPU and GPU of the tasks associated to a specified codelet as a function of the data size of the tasks.

  Codelet copyOverlap West.

Figure 12 .

 12 Figure12. Average of the execution times on CPU and GPU of the tasks associated to a specified codelet as a function of the data size of the tasks.

  eager scheduler.

  b) dmda scheduler.

Figure 13 .

 13 Figure 13. Comparison of the activities obtained with the eager and the dmda schedulers during the computing time.

(a) 0

 0 GPU. Eager scheduler. Execution time = 32.930s. (b) 1 GPU. Eager scheduler. Execution time = 12.465s.

(c) 1

 1 GPU. Source terms forced on GPU. Execution time = 5.810s.

Figure 14 .

 14 Figure 14. Gantt chart of the tasks distribution for different configuration. TOP: no GPU is activated. MIDDLE: one GPU is active and tasks are distributed thanks to the eager scheduler. BOTTOM: one GPU and dmda scheduler.

Table 1 .

 1 Computational time (in seconds) of the integration of system (17) on a 200 × 200 domain divided into 2 × 2, 2 × 4, 4 × 4 or 4 × 8 subdomains, when activating or not a GPU accelerating unit and using the eager or dmda scheduler.

	. Our first order

Acknowledgment

The Hodins team would like to thank the following structures for their support during the 2016 CEMRACS session:

• Groupe Calcul and Fondation Hadamard for its financial support,

• IFPeN for authorizing M. Essadki to be part of the adventure,

• Aix-Marseille and Plafrim (Bordeaux) computational platforms on which most of the results have been performed, • ANR Subsuperjet.