A task-driven implementation of a simple numerical solver for hyperbolic conservation laws

This article describes the implementation of an all-in-one numerical procedure within the runtime StarPU. In order to limit the complexity of the method, for the sake of clarity of the presentation of the non-classical task-driven programming environnement, we have limited the numerics to first order in space and time. Results show that the task distribution is efficient if the tasks are numerous and individually large enough so that the task heap can be saturated by tasks which computational time covers the task management overhead. Next, we also see that even though they are mostly faster on graphic cards, not all the tasks are suitable for GPUs, which brings forward the importance of the task scheduler. Finally, we look at a more realistic system of conservation laws with an expensive source term, what allows us to conclude and open on future works involving higher local arithmetic intensity, by increasing the order of the numerical method or by enriching the model (increased number of parameters and therefore equations).


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 2. 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 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 ) 2 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)

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 (x, y, t n )dxdy. (3) By integrating system (1) over a cell C i,j , one obtains: 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 [10]: whereF (W k,l , W m,n , n) is a chosen numerical flux. In all our computations, the first order local Lax-Friedrichs flux will be used:F Since the local parameter σ LR is equal to the local larger absolute eigenvalue, the numerical scheme (4) is stable in the L ∞ norm, 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)/2.
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). Variable time stepping brings additional problems that are not technically insurmountable, but that complicates the implementation and its analysis a lot. In order to keep it simple and focus on the problems of our interest (task execution, scheduling, etc), we have preferred to fix the time step during the whole simulation.

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 [7]. 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: 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 computationally 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. Now, we detail the glossary attributed to such a new implementation framework within the context of the StarPU runtime.

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 expected 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 resources, following the available kernel implementations. Note that building the dependency graph and launching the task can be an intricate 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 scalability. 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 scalability 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.1). 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.
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 × nVar (ovlpS and ovlpN) and two vectors of size NyLoc × nVar (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 negligible 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.

East West
Copy Copy 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.
• 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. • 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: • 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 (7). 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.

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 to 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 to each variable of each cell and we add RHS to u, see (7). The memory accesses are also coalescent. • checkTimeStep: we associate a virtual processor to 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 and we check if (∆t) i,j is greater than the constant time step ∆t. • internalResidual and borderResidual: we associate a virtual processor to 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 virtual 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 performances.

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 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 output 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. 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 dependency 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 nVar); 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.

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/ρ.

Validation test cases
On this set of equations, we use two classical test cases, namely a linear advection of a density perturbation and a homoentropic vortex. Both have the advantage to present an analytical solution.
At any computational time t, the exact solution is given by 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 homoentropic and its velocity field is a superposition of a constant and a divergence-free fields: From these two last statements, it is obvious that the entropy and the density are simply advected: 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: This is verified on R 2 by a homoentropic rotating Gaussian field: where r is the same radius function as in (10) and R is a characteristic radius of the central perturbation andω is the vortex intensity. Therefore, the initial condition is given by u 0 =ũ +ū, ρ 0 (r) =ρ(r) P 0 (r) =P (r).
with r defined as in (10). The quasi-analytical solution is given by (15) with the definition of r given by (11).
The numerical values of our setup are: 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 (11) 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, [12,13].

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 on CPUs
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.

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.  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. Note the difference in the memory access patterns between ovlpS/ ovlpNand ovlpE/ovlpW, see section 4.3.
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 11-12, 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 11-12, 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 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.   6. 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, [4]: 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. The system is closed through a reconstruction of the size distribution which maximizes the following Shannon entropy: The reconstruction of the size distribution allows to close the two terms: m −1/2 = 1 0 S −1/2 n(S)dS and −n(t, x, S = 0).
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 evaporation 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

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 [4] for more details. This allows to estimate the unknown quantities n(t, x, S = 0) and m −1/2 . Our first order 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 costly, especially when compared to the simplicity of the other terms. Indeed, the maximization of entropy (18) leads to the following non-linear system: This system is solved thanks to Newton-Raphson iterations, as is proposed in [11]. This iterative method has a large arithmetic intensity, which explains its computational cost and justifies its execution on GPUs. Next, computational times are compared, when activating the accelerating device or not, and when using one scheduler (eager) or an other (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 [4]. The initial size distribution of the spray correspondent to a rectangular distribution and the spatially distributed is given by ρ 0 (x, y), which is defined in (10). The initial moments of these initial conditions are scheduler part. 0 GPU 1 GPU where (S min , S max ) = (0.25, 0.75).
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 environment 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 heterogeneous 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 evaporating 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.

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, • Aix-Marseille and Plafrim (Bordeaux) computational platforms on which most of the results have been performed, • IFPeN for authorizing M. Essadki to be part of the adventure, • ANR Subsuperjet, for M. Pelletier PhD grant.