Jean Marie

Couteyen Carpaye

Jean Roman

Pierre Brenner

Design and Analysis of a Task-based Parallelization over a Runtime System of an Explicit Finite-Volume CFD Code with Adaptive Time Stepping

Keywords: HPC, CFD, runtime, task-based

FLUSEPA (Registered trademark in France No. 134009261) is an advanced simulation tool which performs a large panel of aerodynamic studies. It is the unstructured finite-volume solver developed by Airbus Safran Launchers company to calculate compressible, multidimensional, unsteady, viscous and reactive flows around bodies in relative motion. The time integration in FLUSEPA is done using an explicit temporal adaptive method. The current production version of the code is based on MPI and OpenMP. This implementation leads to important synchronizations that must be reduced. To tackle this problem, we present the study of a task-based parallelization of the aerodynamic solver of FLUSEPA using the runtime system StarPU and combining up to three levels of parallelism. We validate our solution by the simulation (using a finite-volume mesh with 80 million cells) of a take-off blast wave propagation for Ariane 5 launcher.

Introduction

For industrial applications of numerical simulation, the most common architecture is nowadays clusters composed of SMP nodes of multicore processors. To develop parallel codes on those machines, a common way is to rely on MPI [START_REF] Forum | MPI: A Message-Passing Interface Standard[END_REF]. While it is possible to use only one core per process and to rely only on Flat-MPI, this approach does not generally scale and it is even worse for codes with a large potential imbalance during execution. A common way to reduce the number of processes while using the same number of cores consists in using OpenMP [START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF] inside the SMP nodes, leading to a two level parallelism.

Another problem is the increasing of the heterogeneity of architectures. Accelerators (e.g. GPGPU, Xeon Phi) are now available and the design of efficient industrial applications exploiting distributed heterogeneous systems with non uniform memory accesses is a complex challenge at large scale. Therefore, applications tend to evolve slowly compared to architectures and achieving performance with a new architecture is a time-consuming task for the developers.

Task-based programming is a good candidate to deal with those issues: describing the problem in a generic manner as a DAG of tasks allows more potential flexibility to exploit the architectures. A runtime system is then in charge of mapping tasks among computational resources (CPU cores and/or accelerators) and of managing memory transfers. By using a powerful abstraction of parallel codes and an efficient runtime system, one can expect to achieve performance quickly on different kinds of architectures (performance portability issue [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures, Concurrency and Computation: Practice and Experience[END_REF], [START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF]).

The aerospace industry faces a lot of problems involving complex unsteady computations for which actual experiences are not doable: e.g. take-off blast wave propagation, stage separation, stability at reentry into Earth atmosphere. Those problems are time-dependent and involve strong shocks.

In this paper, we describe the "taskification" using the runtime system StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures, Concurrency and Computation: Practice and Experience[END_REF] of the aerodynamic solver of the FLUSEPA code [START_REF] Brenner | Three-dimensional aerodynamics with moving bodies applied to solid propellant[END_REF] which is an MPMD MPI/OpenMP code. This aerodynamic solver uses a finite-volume (FV) discretization and an explicit temporal adaptive time stepping scheme. This kind of scheme is well suited for our class of problem because it is conservative and consistent in time ( [START_REF] Kleb | Temporal adaptive Euler/Navier-Stokes algorithm involving unstructured dynamic meshes[END_REF][START_REF] Löhner | Finite element methods for high speed flows[END_REF][START_REF] Pervaiz | Temporal and spatial adaptive algorithm for reacting flows[END_REF]). The method is designed to minimize the computational cost, but it introduces several difficulties for an efficient parallelization (synchronizations, load balancing problems). In this study, we show how a task-based solution over a runtime system can be efficient to tackle these problems.

The paper is organized as follows. Section 2 briefly presents the main computational methods used in FLUSEPA and the existing MPI/OpenMP code. Section 3 presents the runtime system StarPU and Section 4 focuses on the task parallelization of the aerodynamic solver. Section 5 presents an experimental study from an industrial test case, the take-off blast wave propagation. Finally, Section 6 gives some conclusions and perspectives for this work.

Main numerical methods used in FLUSEPA

In this Section, we present the main numerical methods involved in the general computation performed by the FLUSEPA code with a more detailed algorithmic description of the aerodynamic solver on which we will focus in the rest of the paper.

Finite-volume aerodynamic solver

The finite-volume method [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF] is the discretization technique used by FLUSEPA. It is mainly an integral formulation of the conservation laws which are discretized directly in the physical space. In FLUSEPA, we look for the numerical solution of the compressible Navier-Stokes equations in Reynolds-Averaged form (RANS equation for 3D compressible unsteady and reactive flows involving bodies in relative motion [START_REF] Germano | From RANS to DNS: towards a bridging model[END_REF]) that can be written as :

d ∂t ΩCV wdΩ = - ‹ A CV FndS + ΩCV SdΩ (1)
where Ω CV is a fixed control volume (3D cell) with boundary A CV (2D faces), n is the outer-oriented unit normal, w is the conservative variable vector, F is the flux density and S is the source term vector. So, the solver mainly manipulates cells and faces. Field values (e.g. pressure, temperature) are computed for cells and flows are evaluated between faces of cells (refer to [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF] for more details). The main advantage of this discretization method is its conservativity leading to a very close approximation to the physics ; moreover, by using an upwind methodology (Godunov's method [START_REF] Godounov | Résolution numérique des problemes multidimensionnels de la dynamique des gaz[END_REF]) for solving the Riemann problem, one can compute accurately numerical fluxes even when strong shock waves are involved.

To take into account the mesh motion during stage separation, an ALE formulation is used [START_REF] Hirt | An arbitrary lagrangian-eulerian computing method for all flow speeds[END_REF] and when computing fluxes, the intrinsic velocity of each cell face is taken into account. The computation mesh is obtained by a geometric intersection of several meshes (one for each body in relative motion). This conservative technique does not make use of numerical interpolation unlike standard CHIMERA methods [START_REF] Gillyboeuf | Two new Chimera methods: application to missile separation[END_REF][START_REF] Kao | Grid adaptation using Chimera composite overlapping meshes[END_REF]. The volume displacement is also integrated in the associated computational scheme and this is performed by applying the Geometric Conservation Law (GCL) [START_REF] Brenner | Three-dimensional aerodynamics with moving bodies applied to solid propellant[END_REF]. This conservation law is useful in order to transform the geometrical intersection problem into an evolution one. When the relative motion between two meshes is "slow", computing the intersection is not mandatory : the evolution of the cell volumes can be evaluated precisely with the desired level of accuracy through the use of the GCL. Figure 1 shows a separation stage computation: each booster and the main stage are meshed separately. The two main computation steps performed in FLUSEPA at each general iteration, the aerodynamic solver and the mesh intersection computations, are linked by a kinematic computation. During aerodynamic computations, external forces apply to the bodies. From the computation of these loads, a kinematic is obtained. If the bodies moved sufficiently since the last mesh intersection computation, the meshes are moved according to the computed kinematic. Otherwise, the GCL is used in order to take into account the displacement without computing a full intersection between meshes and the kinematics is conserved. This general computational framework is summed up in Algorithm 1. Body displacement computation using the GCL 8: end if

Algorithm 1 General iteration

Temporal adaptive explicit solver

The aerodynamic solver used is explicit and based on a temporal adaptive time stepping scheme. When using an explicit temporal formulation, the maximum allowable time step of a cell is given by its CFL number which depends mostly of the volume of the cell. For an explicit solver, the CFL must be inferior to 1. In classical explicit solvers, the time step is determined by the slowest cell (the cell which has the slowest time step while respecting the CFL condition).

Because we consider complex real size problems, the mesh resolution is not uniform, so using the lowest physical time step would be very penalizing for larger cells. The temporal adaptive algorithm allows to compute each cell near its maximum allowable explicit time step, ranking them in levels; these temporal levels τ are numbered from 0 to a given value θ.

The temporal adaptive method is described in Algorithm 2. At line 1, the maximum allowable time step is computed for each cell. The slowest cell is also found and defines ∆t, the minimum time step in the computation. At line 2, according to ∆t and its maximum allowable time step, each cell is classified inside a temporal level. Inside a temporal level τ, cells are computed at the same time step which is 2 τ * ∆t.

An iteration of the algorithm is composed of multiple subiterations. There are 2 θ subiterations, θ being the level of the fastest cells which are the ones that need only one subiteration to reach the time of the end of the iteration.

The levels τ that are computed are determined by lines 5 to 9. For example, with θ = 3, τ will take successively the following values : 3, 0, 1, 0, 2, 0, 1, 0. Figure 2 shows how level τ evolves after each subiteration in this case; the Predictor: for τ = τ to 0 do

25:

End of predictor for τ : To stay consistent in time, the computations have to be performed in a specific order. When computing a flow between two cells, they must be at the same time. The cells can only have neighbor cells of the same level τ or with the level values τ -1 or τ + 1. If cells are near other cells with a different temporal level, they are positioned at a time that will ensure a consistent computation (Lines 18, 22, 30, 34) [START_REF] Kleb | Temporal adaptive Euler/Navier-Stokes algorithm involving unstructured dynamic meshes[END_REF]. In this way, the computation order is strict and each temporal level is integrated at a specific moment.

The interest of the method depends strongly of the distribution in temporal levels. Let us denote by Ω(τ) the set of cells at temporal level τ and by Ω all the cells in the global domain.

The computational cost of a temporal level τ can be estimated by C(τ) = 2 θ-τ * |Ω(τ)| and the cost of an iteration described by Algorithm 2 is θ τ=0 C(τ). If we compare the computational cost needed to reach the same time with a global time step 2 θ * |Ω|, the cost ratio is

2 θ * |Ω| θ τ=0 2 θ-τ * |Ω(τ)| which is superior to 1.
However, this estimation is an upper bound because of the overhead induced by the temporal adaptive method. Several interpolations are not taken into account in this cost ratio, while their importance increases with the number of temporal levels. More temporal levels imply more cells concerned by interpolation and the cells with a small time step are interpolated more often.

MPI-OpenMP version of FLUSEPA

FLUSEPA uses several kinds of processes to handle the numerical coupling between aerodynamic computations and body movements. Three kinds of processes are used: one process is in charge of coordinating the simulation, some other processes are used to compute the aerodynamic solution and the last kind of processes is dedicated to compute mesh intersections [START_REF] Brenner | Simulation d'interactions aérodynamiques instationnaires autour de plusieurs corps en mouvement relatif[END_REF]. In the following of the paper, we focus on the aerodynamic solver computations of FLUSEPA simulations.

The code uses unstructured meshes in order to take into account complex geometries for interstages. The elements manipulated by the solver are the cells and the faces between them. The current parallel version is based on a two-level parallelism: MPI [START_REF] Forum | MPI: A Message-Passing Interface Standard[END_REF] processes associated with a spatial decomposition and OpenMP [START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF] parallelism inside them.

Making a MPI version of a FV code is usually done by a Domain Decomposition approach. Each process is in charge of a portion of the initial domain and ghost cells are used in order to ensure efficiently communications between subdomains. At the border of a subdomain, faces are duplicated, but each cell belongs to only one subdomain. Figure 6 illustrates this spatial decomposition with 2 subdomains: for the red subdomain, the light red part corresponds to the border cells, the dark red part to inner cells, the light green cells are the ghost cells of this red subdomain, and finally purple faces are the duplicated border faces. MPI communications are only done when necessary according to the temporal level of cells that is currently computed. Border cells are tagged and are computed as soon as possible according to their temporal level. MPI asynchronous communications are used to ensure a good computation-communication overlapping.

The second level of parallelism is achieved in shared memory by using OpenMP directives (OMP DO) applied to loops concerning the cells and faces inside a subdomain.

The main problem of this spatial domain decomposition is the fact that the cells have not the same computational cost which is determined by their temporal level. To ensure a better load balancing, we give a weight to each cell of the subdomain according to its temporal level. However, the temporal level of a cell can change between iterations and a recomputation of a new domain decomposition is needed periodically. Despite that, the way the time is integrated leads to an important time wasted in synchronizations. The time integration implies a strict order for the cells to be processed depending on their temporal level: neighbor cells must be at the same time during computation. This temporal locality information is partially lost with the current parallelization.

Figure 3 shows an execution trace of the MPI version using 8 processes for θ = 4 (16 subiterations). We performed an instrumentation of the code allowing to identify each subiteration; in this trace each one is colored with a different color, the first subiteration being colored in green. When a subiteration concerns only τ 0 cells (one over two), it is colored in cyan. Time spent in MPI functions is identified in red. Computation not concerning the solver is colored in black. After the last operation in black, the process has no more work to accomplish for the current iteration. Process 4 conditions the time of the iteration, being the one that finishes last. Most of the cells of a temporal level inferior to 4 are computed by this process 4. At the first subiteration, the only one during which τ 4 cells are computed, a lot of time is wasted in synchronization. For other subiterations, this process 4 computes without any interruption. Another interesting point is that processes 5, 7 and 8 can start the 9th subiteration (in yellow) before the others. Cells processed during this subiteration have τ 3 ,τ 2 , τ 1 and τ 0 temporal levels. This means that τ 3 cells of those processes are not related to the one of the other domains (they do not share borders). Other processes (1, 2, 3, 4) can compute this subiteration lately.

This synchronization issue is one of the key elements to justify the development of a task-based version of the aerodynamic solver. By working on subdomains inside each process, we want to be able to capture all the dependencies during computations and to exploit more asynchronism with the help of a runtime system.

The StarPU task-based runtime system

There exist different libraries and frameworks to exploit task-based parallelism: e.g. SMPSs [START_REF] Duran | Extending the OpenMP tasking model to allow dependent tasks[END_REF], StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures, Concurrency and Computation: Practice and Experience[END_REF] , PaRSEC [START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF], CnC [START_REF] Budimlić | Concurrent collections[END_REF], Legion [START_REF] Bauer | Legion: Programming distributed heterogeneous architectures with logical regions[END_REF], SuperMatrix [START_REF] Chan | Supermatrix out-of-order scheduling of matrix operations for SMP and multi-core architectures[END_REF].

As said in the introduction, the main goal of task-based programming is to ensure performance portability on heterogeneous manycore distributed platforms. In this framework, the algorithm is described as a sequential task flow with data dependencies expressed through read/write attributes for each task parameter. This task flow is translated into a Direct Acyclic Graph (DAG) of tasks: the nodes represent the tasks and the edges between the nodes are the dependencies. Then, a runtime system is in charge of scheduling the tasks over the computational resources (CPU, GPU) and of managing the data transfers.

Good results have been achieved by this approach for dense linear algebra [START_REF] Bosilca | Flexible development of dense linear algebra algorithms on massively parallel architectures with DPLASMA[END_REF] and sparse linear algebra [START_REF] Lacoste | Taking advantage of hybrid systems for sparse direct solvers via task-based runtimes[END_REF]. There are some other works about parallelization of applications over runtimes: S3D over Legion [START_REF] Bauer | Legion: Programming distributed heterogeneous architectures with logical regions[END_REF], ScalFMM over StarPU [START_REF] Agullo | Task-based FMM for multicore architectures[END_REF].

To construct the DAG, there are two common approaches. First, a Parametrized Task Graph (PTG) can be used [START_REF] Cosnard | Automatic task graph generation techniques[END_REF]: tasks are not enumerated but parametrized and dependencies between tasks are explicit. Another way is to use the Sequential Task Flow (STF) model. With STF, tasks are inserted from the main program and the dependencies are computed at task insertion according to data accesses [START_REF] Allen | Optimizing compilers for modern architectures: a dependence-based approach[END_REF]. One advantage of the later model is the fact that tasks can be inserted according to the results of previous computations.

The StarPU runtime system [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures, Concurrency and Computation: Practice and Experience[END_REF] relies on the STF model and computation resources (e.g. CPU, GPU) are seen as workers.

Tasks are inserted from the main program through calls to starpu_insert_task. The dependencies between the tasks are then computed by the runtime system according to the data accesses and how they are accessed (read, write, read-write). This implies that the DAG is unrolled during the execution. The task insertion is asynchronous: the computation can start even if not all the tasks have been inserted. When a task is ready (i.e. all its dependencies have been fulfilled), it becomes available to the scheduler. Then, according to the scheduling strategy, workers select tasks and execute them. Some hints are available to the scheduler: at task insertion, it is possible to give a priority to a task and some schedulers can take advantage of this information. It is also possible to have performance models which will compute a weight for a task according to various parameters (e.g. size of the data, targeted architectures) in order to help the scheduler. In StarPU, the scheduler is a modular component : it is possible to define multiple scheduling strategies. In order to do so, the programmer has to define containers for the tasks, a push operation that will put ready tasks inside the right container and a pop operation which will allow workers to choose a task to be executed. Figure 4 shows how the built-in scheduler "prio" works. Circles represent tasks, the number inside represents the priority associated. Tasks are inserted with a priority from 0 to 4. When they are ready, they are pushed to the right container. The containers are FIFO lists : one list exists for each allowed priority. For the pop operation, the workers will check every list by order of decreasing priority until they find a task to execute.

Parallel tasks and worker-contexts

It is possible to exploit existing OpenMP code within StarPU through the use of context [START_REF] Cojean | Exploiting two-level parallelism by aggregating computing resources in task-based applications over accelerator-based machines[END_REF]. A context can be seen as a set of computational resources with a scheduling strategy in which tasks can be submitted. Contexts can be nested.

When no scheduling strategy is specified, the context is seen as a worker (called worker-context): when a task is pushed on this worker-context, it is supposed to be executed on all its computational resources.

To use an OpenMP code inside a worker-context, a specific initialization task that will bind OpenMP threads to the CPUs of the worker-context must be used. By default, tasks are inserted in a global context. So, if one wants to create 4 worker-contexts of 4 CPUs, he has to create a main-context with a scheduling strategy that will contain these worker-contexts and insert tasks in this main-context. There are two ways to use MPI parallelism with StarPU. First, StarPU provides starpu_mpi_insert_task that lets StarPU handle all the communications and data transfers. With this approach, all the tasks must be inserted in all the nodes of the cluster. Communications are implicit and each node knows which node will execute a particular task.

Distributed parallelism with StarPU

starpu_mpi_isend_detached and starpu_mpi_irecv_detached allow to describe communications explicitly. When using those primitives, communications are consistent with the computation dependencies inferred by task insertion inside the nodes.

Figure 5 shows the difference between the two approaches while using three nodes. The same DAG is considered. In the implicit case, the whole DAG is present in each node and only a portion of it is useful. In the explicit case, each node only contains the tasks it will execute and communication tasks are present.

MPI + Task design of the aerodynamic solver

Computation Elements

As the first level of parallelism we can exploit in FLUSEPA is spatial, the natural way to described tasks is to use multiple subdomains. In the DAG, we want to express dependencies in such a way that allows the maximum concurrency. In order to achieve this, we use an algorithmic abstraction called "Computation Element" (CE).

The code mainly manipulates faces and cells and this must be exploited to have a good task parallelization. Four main basic computation patterns can be found in the original code: computation on cells, computation on faces, computation on faces using cell values, computation on cells using face values. However, the way the code has been written does not lead to a clear distinction between those patterns when they were used, and thus a rewritting step has been performed. Most of the computation kernels of the aerodynamic solver will now correspond to one of these computation patterns in order to achieve a well-structured task version (cf. Algorithm 3).

The domain partitioning is done with the graph partitioner SCOTCH1 [START_REF] Pellegrini | Scotch: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs[END_REF]. For each subdomain created, we associate a CE containing informations to retrieve each field or flow value (See Section 2.1). The different components end for 12: end for 13: Wait all tasks 14: Send informations to master process 15: Foreach cells update intensive values 16: Wait all tasks In order to generate tasks, we mainly use the different data types we described previously and exploit the associated computation patterns. We rely on the STF model of StarPU. For each pattern, we used a "foreach" function that is in charge of generating the right tasks. We mainly exploit Algorithm 2, each line being converted in one or several "foreach".

The task insertion follows Algorithm 3. The time step computation (line 1) and the classification of cells in temporal levels for each CE (line 2) are also done using tasks. Currently, we have one synchronization after the classification of cells in temporal levels. In order to know which task must be inserted, the temporal levels inside each CE must be known. Some CEs do not contain cells of some temporal levels and this fact is exploited to avoid insertion of useless empty tasks.

The foreach functions can generate a high number of tasks and in particular chains (succession of tasks that depend only of one previous task). In order to reduce the number of tasks, we implemented a strategy to pack tasks at runtime. A pack is then a single task composed of multiple elementary tasks.

Computation kernels are written in such a way that they only write one component of a CE and the pack mechanism relies on that. We have 3 kinds of task packs, each one for modifying data either on border cells, inner cells or faces. For each CE, these task packs exist.

The foreach functions know which type of component (face or cell) will be modified. If it is the same at the previous one, tasks will be added to the corresponding pack associated to the CE. Otherwise, the previous packs are inserted as tasks and new packs are created.

Another case that generates chains of tasks happens when cells of a given temporal level are present only in the inner cell component of a CE. For each CE, we check if it is the case at line 4 of Algorithm 3. To handle this particular context, we introduce a 4th kind of pack called "large task pack". As there is no interaction with other CEs, they contains all the tasks for a given temporal level and a given CE, disregarding the type of component currently modified.

Optimizing this point is critical because tasks that work on cells with low temporal levels are numerous and they work on a low number of cells. This fact is shown in Table 1 which gives the proportion of cells in each temporal level and the proportion of associated computation cost for the test case we use in Section 5 when fixing the maximum temporal level to 2. The large majority of cells are of the higher temporal level (τ 2 ). Few cells belong to low temporal levels, so we can expect that most of them won't be part of the border component of a CE.

τ 0 τ 1 τ 2
Cells 0.05% 2.42% 97.53% Computation 0.20% 4.72% 95.08% The colored diamonds correspond to computations that modify data for faces of a CE, black diamonds for inter-CE faces, small circles for border cells and large circles for inner cells. Triangles represent large task packs: in this case, each large task pack contains more than 10 subtasks. In this DAG, we can notice that τ 0 temporal level cells are only present in the inner component of the red CE as there are large task packs. Increasing the number of CEs would enlarge the width of the DAG and using more temporal levels would enlarge its height for some CEs. The DAG can really be unbalanced according to the cell distribution in temporal levels inside each CE, so the way the graph is traversed can have a strong impact on the efficiency of the computation.

MPI + Task parallelization

As we insert tasks according to previous computations, we decided to rely on explicit communications. Each process gets a domain after the domain partitioning. At this moment, for each process, border cells are identified and border faces are duplicated. In order to distinguish border cells of the domain associated to a process and border cells of the component of a CE, we note the first ones "MPI-border cells".

These MPI-border cells are then included to the border component of the CE they belong. After that, each process communicates with its neighbors to complete informations for the common border cells. For each (local CE, foreign CE) couple, we create a border ghost cell component.

Figure 8 shows two domains with two CEs each. For the two CEs that share a border, we can see the new components : MPI-border cells, inter-domain faces and ghost cells.

For the management of communications, we use for sending data a task that copies border cell content to a temporary buffer and this buffer is sent using starpu_isend_detached. For receiving data, the starpu_irecv_detached operation is used in a temporary buffer and then, we copy the data in the ghost cell component.

We use a temporary buffer because we do not want to send the whole border component of a CE: some border cells are not MPI-border cells and when processing low temporal levels, a fraction of the MPI-border cells is concerned.

Before foreach functions that uses cells in read mode, we insert communications as tasks.

Our task-based implementation is able to exploit 3 levels of parallelism : between computation nodes, MPI is used; inside a SMP node, StarPU tasks are used; and tasks can be parallel by using OpenMP.

Experimental study

For the experimental validation of this work, we used a cluster with nodes composed of two 8-cores Sandy Bridge and 64 GB of RAM.

The test case from Airbus Defence & Space is the computation of a take-off blast wave. For an Ariane 5 rocket, boosters provide 90% of the thrust at lift-off and the objective is to compute the resulting blast-wave implied by the ignition of the boosters. The mesh is composed of 10M cells and is finer around area of interest (Figure 9). During take-off blast wave propagation, there are two overpressures: the first one is due to ignition of the boosters, the second is due to ducts. Those events are visible at Figure 10: the comeback of the wave can be seen from the 4th picture.

We consider a case with θ = 4. In this computation, the cell repartition in different temporal levels and the theoretical computation cost associated evolve during the computation, which is shown at Figure 11 and Tables 2 and3. The tables show the proportion of cells and the associated computation cost at two different times of the computation: Table 2 at the start of the computation and Table 3 at the 2300th iteration.

Figure 11a shows how the proportion of cells in each temporal level evolves. There are few cells of the lower temporal level τ 0 and there are not visible on the figure. We can observe that the proportion of cells of temporal level τ 1 evolves during the computation, but they represent at most 2% of the cells. Figure 11b shows the relative cost of each level τ and its evolution during the computation. Figure 11c shows the global cost of an iteration (in floatingpoint operations) during the computation. It shows that the cost of one iteration varies and that it gets more and more expensive while the computation progresses. The computation cost associated to τ 1 cells evolves a lot. According to Tables 2 and3, it increases from 0.45% to 12.45% between the start of the computation and the 2300th iteration.

For our task-based parallel implementation, we focus now on one iteration of the aerodynamic solver and we first study in shared memory the impact of several parameters: the priority strategy for the scheduling of tasks, the number of CEs, the number of parallel workers. Once we have evaluated the impact of those parameters, we will evaluate the implementation at two different iterations with different cell proportions in the temporal levels. Then, we will evaluate a distributed memory version combining all the levels of parallelism. 3: Cell distribution and associated computation cost at the 2300th iteration of the computation (θ = 4).

Shared memory study

Impact of priorities and of number of CEs

For this study, we don't use parallel workers, so each CPU is associated to only one worker.

We work with the StarPU built-in "prio" scheduler. In this scheduler, tasks are pushed in different priority queues and they are popped by priority order.

We consider two strategies. In the first one, we don't give any priority to any CE. For the second one, we give to tasks a priority according to the CE they belong. We want to give an advantage to the CEs with low temporal levels. However, when giving a priority to a task in StarPU, it doesn't automatically propagate the priority to the predecessors of the task in the DAG. So, it is also necessary to give a good priority to the predecessors of the tasks.

To achieve this goal, we give the highest priority to CEs which contain cells of temporal levels 0 or 1. Then for each CE, we compute a priority : we evaluate the lower distance from this CE to the other CEs previously prioritized. Then we map priorities according to the distance computed: the lower the distance, the higher the priority.

The result for our test case are shown at Figure 12. The reference time (20.02s) is the one for a computation with only 1 CE and one unique parallel worker which uses the 16 cores of the node; this is the configuration that mimics the best the previous OpenMP version in shared memory. The first observation is that this reference configuration with 1 CE is the worst in terms of performance. For a number of CEs from 16 to 256, we compare the elapsed time with our priority strategy (on the right) and without priority (on the left).

We show the state of the different workers during one iteration. The global size of the bar indicates the time (in seconds) needed to complete one iteration and the fill colors correspond to the proportions spent in each state (executing, sleeping, overhead). The overhead state contains, among other things, time spent to compute dependencies and insert tasks. The sleeping state means that the worker is ready but there is no ready task available.

The TS/GER part corresponds to lines 1 to 3 of Algorithm 3, including the "Wait all tasks" step. The Solver part corresponds to lines 4 to 15. We observe that the time spent in computation (blue and green bar) does not evolve 13 when we increase the number of CEs, or if we use or not the priority strategy. Concerning computations from 16 to 128 CEs, sleeping state (red and dark pink bar) is reduced by decomposing in more CEs and by using the priority strategy. This strategy is beneficial as soon as we use 32 CEs. Concerning the overhead (yellow and gold bar), it evolves linearly with the number of CEs and the priority strategy has no influence. We detail below what happens for 128 and 256 CEs. The first step of the computation is in black (TS/GER, lines 1 to 3 of Algorithm 3). For the solver part of the computation (lines 4 to 13), tasks are colored according to the subiteration they compute.

Of course, the same DAG is built and only the order of computations changes when using the priority strategy. In the trace without priorities, we observe that at the end of the computation, there is a large starvation zone: workers spend time in sleeping state (in red) because there is no task available. When using priorities, this starvation is much less important: we perform the DAG traversal in a way that allows more tasks to be available at the end of the computation. It is also noticeable that we managed to finish the scheduling by task concerning the first subiteration.

The evolution of ready tasks in the solver is shown at Figure 14. For the two strategies, we represent the number of ready tasks over time. Ready tasks are the ones already available for the workers: their dependencies have been fulfilled and they can be executed. The spike in the end comes from line 15 of Algorithm 3: for each cell component of a CE, a last task is inserted, just after a "Wait for all" (line 13). Some operations that are not done using task for legacy reasons are performed line 14. Traversing the graph without priority unlocks more tasks at start, but in the end, not enough tasks are available to feed the workers. Our strategy manages to keep more ready tasks available until the end.

The idle step we observe in the TS/GER part of the computation is due to the fact that with our task granularity, we execute tasks quicker than we insert them. Finally, we can notice that the gain achieved by using the priority strategy is almost 10% (14,18s versus 15,56s). When using now 256 CEs, we notice that the sleeping time of both the solver and the TS/GER parts increases. This can be explained by the time needed to insert tasks. For example, it takes 13.6 s to insert all the tasks of the solver step (from lines 5 to line 12) and the total computing CPU time is 195.8 s (about 12.2 seconds per worker). So, increasing the number of CEs allows more concurrency and the idle time is reduced. But with 256 CEs, the task insertion becomes so costly that it affects the global computation performance.

In conclusion of this first study in shared memory, we can say that the way we prioritize the tasks is efficient: we manage to feed the workers as long as possible whereas starvation was important without priority. The task-based description truly allows to take advantage of the irregularity of the computation.

Study with parallel workers

In order to exploit all the CPUs without creating too much tasks, another way is the use of parallel workers. Instead of being executed on only one CPU, a task is then executed on several ones. We rely on OpenMP DO loops for the parallelization of our computational kernels, which comes from the previous OpenMP version of the code.

We consider different numbers and sizes of workers (the size is the number of cores used), from 16 workers of size 1, to one worker that uses all the cores: so we have the relation "(number_of_workers)×(size_of_workers)=16". CPUs that belong to the same worker are on the same socket, except when the 2 sockets are used by one worker. We test configurations from 16 CEs to 128. The previous strategy for scheduling tasks is always used in this study.

At Figure 15, we first notice that the total time spent in computation tends to increase slightly while we use fewer and larger workers. Indeed, the scalability of our computational kernels is not good enough. However, this is not critical until a worker size equal to 8. When we use the whole node and its two sockets to form only one worker, the situation is more critical. The overhead evolves linearly with the size of workers, while the idle time is reduced when using parallel workers. The best overall configuration is obtained for 32 CEs and 8 workers of size 2.

Different cells repartition in temporal level

At Figure 11, we highlighted the fact that the cell repartion evolves and so the associated computation cost during the simulation. In Tables 2 and5, we detailed two specific iterations, one at the start of the computation (iteration 1) and the other at the 2300th iteration. The main difference between the two cases is the computation cost produced by the cells of low temporal levels. Notably, for τ 1 cells, the computation cost is less than 0.5% in the first iteration and 12.45% in the 2300th.

In order to see the impact of this change, we compare now the behaviour of the implementation at these two iterations. We use 3 different parallel configurations (16 × 1, 8 × 2, 4 × 4) and 4 different numbers of CEs [START_REF] Duran | Extending the OpenMP tasking model to allow dependent tasks[END_REF]32,64,128). The "prio" strategy is always used.

For the first iteration, Table 4 gives some informations about the DAG that have been generated. The Tsub column indicates the time needed to insert all the tasks, the elementary tasks column indicates how many tasks would have been generated without the usage of pack, while the inserted packs column indicates how many tasks were actually inserted inside the runtime. Figure 16 shows the performance obtained. The black line shows the performance for the OpenMP version. Table 5 and Figure 17 give the same informations for the 2300th iteration. If we look at both Tables 4 and5, we can notice first that the packing strategy tends to reduce the number of tasks that the runtime has to manage. Increasing the number of CEs raises the number of tasks and the time needed to submit these tasks. If we compare the number of taks inserted between the two cases, we see that the number of tasks inserted is more important at the 2300th iteration. Considering now Figures 16 and17, it is important to notice that the 2300th iteration is more expensive than the first one. This is true for both OpenMP (25.88s versus 20.32s ) and all task-based versions and this is consistent with Figure 11c.

Again, we can notice that adding more CEs reduces the synchronization time. Using larger workers reduces the synchronization time but has a negative impact for the time spent executing tasks and increases the overhead. But if we look at the synchronization time for the 2300th iteration, it is more important than for at the first one. With 16 CEs, only the 4 × 4 configuration performs better than the OpenMP version. Without parallel workers (configuration 16 × 1), the synchronization time in this case is more important than the time spent executing and so the task-based version is beaten by the OpenMP one. With 32 and 64 CEs, the synchronization time is still important and the taskbased version without parallel workers is beaten by the OpenMP one. This result emphasizes the gain due to the use of parallel workers : the 4 × 4 configuration performs better than all the others from 16 to 64 CEs and is only beaten by the 8 × 2 one with 128 CEs. As shown by Tables 4 and5, increasing the number of CEs implies a more expensive task submission cost. The impact of parallel workers on this submission cost is neglectable and their use allows to greatly reduce the time spent in synchronization.

A distributed memory version

Using the task description, we hope to reduce the synchronization time we noticed with the initial MPI+OpenMP version and that is higlighted by Figure 3. We have shown that the task based version trully allowed to compute in a more relaxed order. This is notably shown in Figure 13 where we can see that tasks corresponding to the first subiteration can be executed at the end of the iteration.

To evaluate the task version in distributed memory, we first start with a study using the same test case than previously in shared memory. Then, we use another mesh with a simpler geometry but with a greater number of cells. In each case, we will compare the task version with the initial MPI+OpenMP version.

First test case

We still use the same test case as in the shared memory study. We consider here experiments from 2 to 8 MPI processes, with 8 to 32 CEs per process, and we configure our workers in 3 configurations : 8×2, 4×4 and 2×8. Each process is in charge of a portion of the mesh and the domain decomposition takes into account the cost of cells. Figure 18 shows the different elapsed times for the different configurations tested. We can see that the best performance (4,331s) is obtained with 8 CEs per process and 6 processes configured as 4×4 (4 workers and 4 cores per worker). This case is detailed in Table 6. The load balancing is good with a difference of 13% of executing time between Process 1 and Process 4. However, the idle time is important: with only 2 CEs per worker, it is not easy to exploit enough asynchronism and communication-computation overlapping.

Figure 18d shows the elapsed time for the initial MPI+OpenMP version of FLUSEPA. We consider 2 configurations, one with one process per node, and one with one process per socket. Unfortunately, the problem does not fit in memory for a larger number of processes per node, mainly because the memory consumed by process depends of the load balancing. We see that the absolute best performance is obtained for this previous version with 2 processes per node and 10 nodes. But when we consider a number of nodes for which the task-based version is competitive (at most 6), the task-based version gives a better result than the previous one (4.33s versus 5.98s) for 2 processes per node and only 6 nodes.

With 8 MPI processes, we don't reduce the elapsed time in the solver. This is due to the number of CEs which varies from 64 for 8 CEs per process to 256 for 32 CEs per process. When we work on more nodes, the computation time by node decreases, but the overhead stays almost the same. This is the main reason why we cannot use an important number of CEs. Load balancing is also a concern. The size of the problem is a little bit tiny to overcome the overhead induced by the parallelization method.

A larger test-case

In this section, we evaluate the behaviour of the task-based version on a larger problem and with a larger number of nodes; we compare its performances to the ones of the initial MPI+OpenMP version.

The test case is still a take-off blast wave propagation with a simpler geometry but with 80 M cells. Due to the simpler mesh, only 4 temporal levels are present at most in this computation (θ = 3).

The cluster used was composed of 20 cores by node (2 Ivybridges with 10 cores each).

We run the initial MPI+OpenMP version with one process (named MPI-OMP-1) or two processes (named MPI-OMP-2) by node and we use a 4 × 5 configuration for the task-based version (4 workers each using 5 cores). This task-based version is performed by using 10 CEs by process (named RT-10CEs) or 20 CEs by process (named RT-20CEs).

In order to be able to run the computation without memory swapping, we must use at least 16 nodes2 . Concerning the MPI+OpenMP version using 16 nodes, only the computation with one process by node could be achieved due to a memory consumption problem when using 2 processes by node. The computation is done for 16, 20, 24 and 28 nodes and Figure 19 shows the results achieved.

Regarding the MPI+OpenMP version, the configuration using 2 processes by node is clearly better than the one with only one process by node.

The task-based versions using StarPU are clearly more efficient than the MPI+OpenMP versions. With 16 nodes, the best task-based configuration (20 CEs by node) is 70% faster than the best MPI+OpenMP version (with 1 process by node). With 28 nodes, this gain is 49%, but here the best task-based configuration uses 10 CEs while the best MPI+OpenMP version uses 2 processes by node. We notice that the task-based version using 10 CEs by node scales better than when using 20 CEs by node. Indeed, the computation granularity becomes more and more problematic as the number of nodes increases (the number of cells by CE is lower when using 20 CEs by node).

Figure 20 gives some details about the time spent in the different states for the different configurations and for the different numbers of nodes.

The Y-axis shows the total time for an iteration while the bars are colored according to the proportion spent in each state: in red for the time spent in sleeping, in blue for the time spent in computing and in yellow for the overhead.

For the MPI+OpenMP version, we notice that using 2 processes by node reduces significantly the time spent in computing. The proportion of time spent in the sleeping state remains around 50% in all the cases.

If we compare the two configurations used in the task-based version, the best one for 16 nodes is the one using 20 CEs by node. For this number of nodes, increasing the number of CEs allows to reduce the time spent in the sleeping state. However, the overhead is more important than with 10 CEs by node. If now we use more nodes, the best version becomes the one using 10 CEs by node. The overhead implied by the increasing number of CEs become more important. Furthermore, by adding more nodes, we increase the total number of CEs and each CE contains less and less cells.

It would be nice to be able to limit this overhead in order to obtain a better gain when using the task-based version. The easier way to get more asynchronism is to increase the number of CEs, but this also increases the associated overhead which currently limits the potential gain.

Conclusion and perspectives

In this paper, we have described a task-based distributed version of an aerodynamic solver with time adaptive time step. We have validated our implementation with several real-life industrial cases.

When under the right conditions (i.e. when the overhead induced is not too high), the implementation shows very promising results. The parallelism offered by the task-based paradigm allows a better exploitation of the computational resources. On the largest test cases, the runtime version performs much better than the initial MPI+OpenMP version in every situation evaluated. However, it is not possible to remove all synchronizations : if we use too much CEs per process, the overhead induced is too high. Some code analysis should be performed in order to see if it is possible to reduce this overhead.

It is clear that exploiting parallel tasks is useful reducing the number of workers, but this point can be limited by the scalability of the computational kernels. Concerning the way we exploit parallel tasks, we tested only the simple case of multiple workers of the same size. It is possible to explore more heterogeneous configurations (e.g. 1×8+2*2×4 : one worker of size 8, and two of size 4). It could even be possible to modify the configuration on the fly regarding current status of computations. To overcome the starvation seen at the end of the iteration, one possible strategy is to gather all CPUs into one worker when there is few ready tasks. Context resizing is possible in StarPU and have been described in [START_REF] Hugo | Composing multiple StarPU applications over heterogeneous machines: a supervised approach[END_REF].

Right now, the decomposition into CEs is fixed at the start of the computation, but the temporal level of cells evolves during the computation between iterations. Having the possibility to reshape the CEs could lead to an improvement : it should be possible to generate a more optimized DAG with CEs that take into account the temporal levels of cells. Another use for the CE reshape could be load balancing for the distributed version. Between iterations, it could be possible to exchange CEs between processes which would allow a better load balancing. Furthermore, if we manage to remove the hard synchronization we have at the end of general iterations, we can expect to pipeline iterations while applying an asynchronous load balancing scheme.

We intend also to develop a task-based parallelization of the intersection mechanism of FLUSEPA in order to exploit more asynchronism in the whole application when considering booster separation simulations.

This work is a proof of concept showing that a task-based parallelization over a runtime system can be a good approach on the route towards efficient large scale simulations on forthcoming post-petascale supercomputers. Achieving good performances could be more easy from a task-based description : adjusting the grain of parallelism will be possible and in our case, it can be done by choosing the number of CEs or the configuration used for parallel workers. Overlapping of communication by computation can also be taken to another level. First, the fine description of dependencies allows to remove synchronizations. Also, it is possible to have scheduling strategies that favor communication/computation overlapping by selecting the right tasks that will generate communications as soon as possible. However, the programming effort in order to get a task-based version is important. This may be reduced by the rise of programming standards (OpenMP 4.5 [START_REF]OpenMP application program interface version 4.5[END_REF] now allows to do things similar to STF through the depend clause).

Figure 1 :

 1 Figure 1: Booster separation.

Algorithm 2

 2 Temporal adaptive time stepping scheme in FLUSEPA: one iteration of the aerodynamic solver 1: Time step computation 2: Classification of every cell inside a temporal level 3: Temporal adaptive loop: 4: for subiteration=1 to 2 θ do for tmp = 1 to θ do 7:if (mod(subiteration -1, 2 tmp ) == 0) then 8:

  evolving levels τ are in red. The physical time reached after one iteration with θ = 3 is equivalent to 8 iterations with a global time step.

Figure 2 :

 2 Figure 2: Time reached after each subiteration for each level τ (θ = 3; 8 subiterations).

Figure 3 :

 3 Figure 3: MPI trace, 8 processes (θ = 4; 16 subiterations).

Figure 4 :

 4 Figure 4: Scheduler with priorities ("prio").

Figure 5 :

 5 Figure 5: DAG representation inside nodes with implicit or explicit communications.

Figure 6 : 5 : 8 :

 658 Figure 6: Illustration of two CEs.

Figure 7 :

 7 Figure 7: DAG for a computation with 6 CEs and θ = 2.

Figure 7

 7 Figure 7 shows a DAG for an iteration with a maximum temporal level θ = 2, three temporal levels of cells, 4 sub-iterations and 6 CEs. Colors represent the different CEs. The red CE contains cells of the three levels, two other CEs (blue and purple) contain τ 1 and τ 2 temporal level cells. The remaining three CEs contain only τ 2 temporal level cells.The colored diamonds correspond to computations that modify data for faces of a CE, black diamonds for inter-CE faces, small circles for border cells and large circles for inner cells. Triangles represent large task packs: in this case, each large task pack contains more than 10 subtasks. In this DAG, we can notice that τ 0 temporal level cells are only present in the inner component of the red CE as there are large task packs.

Figure 8 :

 8 Figure 8: Two domains, each with two CEs.

Figure 9 :

 9 Figure 9: Mesh for the take-off blast wave computation.

Figure 10 :Figure 11 :

 1011 Figure 10: Take-off blast wave computation.

Figure 12 :Figure 13 :

 1213 Figure 12: Impact of Priority Strategy with varying number of CEs.

Figure 14 :

 14 Figure 14: Evolution of ready tasks (128 CEs).

Figure 15 :

 15 Figure 15: Time performance with different configurations with parallel workers.

Figure 16 :

 16 Figure 16: Performance (θ = 4; first iteration).

Figure 17 :

 17 Figure 17: Performance (θ = 4; 2300th iteration).

(a) 8 Figure 18 :

 818 Figure 18: Time spent in the solver with different configurations (X-axis gives the number of nodes).

Figure 19 :

 19 Figure 19: Comparaison between MPI+OpenMP versions and task-based versions on a large test case.

Figure 20 :

 20 Figure 20: Time spent in each state for a larger test case with different configurations.

Table 1 :

 1 Cell distribution and associated computation cost (θ = 2).

Table 2 :

 2

		τ 0	τ 1	τ 2	τ 3	τ 4
	Cells	0.02 % 2.00 % 3.04 % 5.01 % 89.93 %
	Computation 0.22 % 12.45 % 9.48 % 7.81 % 70.04 %

Cell distribution and associated computation cost at the start of the computation (θ = 4).

Table

Table 4 :

 4 DAG creation information (first iteration).

	#NCEs	Tsub	# elementary tasks # inserted packs
	16	0.7 ± 0.1 s	23506	2826
	32	1.8 ± 0.1 s	43068	5461
	64	3.4 ± 0.3 s	64955	8696
	128	7.7 ± 0.5 s	118633	16259
	#NCEs	Tsub	# elementary tasks # inserted packs
	16	0.6 ± 0.1 s	34140	4238
	32	2.0 ± 0.1 s	58061	7361
	64	3.4 ± 0.3 s	90383	12006
	128	8.0 ± 0.5 s	157605	21174

Table 5 :

 5 DAG creation information (2300th iteration).

  TS/GER idle TS/GER executing Solver idle Solver executing Submission time Elapsed time

	Process #1	0.338	0.197	1.440	2.356	0.608	4.331
	Process #2	0.357	0.174	1.559	2.241	1.301	4.331
	Process #3	0.351	0.184	1.401	2.395	0.937	4.331
	Process #4	0.293	0.243	1.790	2.005	0.218	4.331
	Process #5	0.310	0.216	1.599	2.205	0.599	4.331
	Process #6	0.344	0.191	1.627	2.169	1.259	4.331

Preprint of an article accepted for publication in Journal of Computational Science. (https://dx.doi.org/10.1016/j.jocs.2017.03.008) Accepted manuscript available on arXiv. (https://arxiv.org/abs/1704.01144v2)

Couteyen Jean Marie / Journal of Computational Science 00 (2017) 1-22

https://gforge.inria.fr/projects/scotch/

Table 6: Average time (in seconds) for each process. We use 6 processes, 8 CEs per process, 4 workers and 4 cores per worker.

The cluster used did not allow to swap.