D A Mbengoue

D Genet

C Lachat

E Martin

M Mogé

V Perrier

F Renac

F Rué

M Ricchiuto

Comparison of algorithm in Aerosol and Aghora for compressible flows

Keywords: Finite elements, parallel performance, compressible flows Finite elements, parallel performance, compressible flows

This article summarizes the work done within the COLARGOL project during CEMRACS 2012. The aim of this project is to compare the implementations of high order finite elements methods for compressible flows that have been developed at ONERA and at INRIA for about one year, within the AGHORA and AEROSOL libraries.

Introduction

On the need for high order in aerodynamics Over the last three decades, most of the simulations for designing aircraft have been reduced to second order finite volumes methods, mainly with RANS (Reynolds Averaged Navier-Stokes) models for turbulence. The usual trick for having an accurate solution with a RANS simulation consists in performing a sequence of computations on more and more refined meshes. However, the RANS approach can be inefficient for simulating some problems that are intrisically time dependent, because it is a time averaged model, and thus stationary. For time dependent simulations, the most accurate method is called DNS (Direct Numerical Simulation) and consists in meshing the domain at the turbulent scale, without any turbulence model. However, due to the required mesh fineness, DNS is often considered only for low Reynolds numbers (since the number of cells increases as Re 9/4 ) and/or for academic configurations, for which efficient finite differences methods can be used. An intermediate solution, between DNS and RANS approach consists in applying a spatial filter to the Navier-Stokes equations. The resulting model is time dependent, but shall be closed for taking into account small scales. This method is called the Large Eddy Simulation (LES) approach. For both DNS and LES, the accuracy yielded by adaptive mesh refinement is harder to implement because the problem is time dependent. It would therefore imply theoretically to perform a mesh adaptation at each time step, which is very costly. It also raises issues regarding dynamic load balancing as far as parallel computing is concerned.

How to derive a high order method

Mathematically speaking, a sufficient condition for having a high order approximation of a function is to find a piecewise polynomial approximation of this function, for example by interpolation or L 2 -projection. Once the type of approximation is chosen, a numerical scheme must be derived from this approximation. The following methods can be considered.

• High order finite volume methods

In finite volume approximations, the solution is considered as piecewise constant. For getting a high order method, a high order polynomial representation can be computed, by interpolating the values on the neighboring cells. For example, in one dimension, a second order polynomial can be computed on one cell by interpolating the values on the left and right cells. This can be generalized to higher dimensions and to higher order polynomial approximations, and the price to pay for having a high order representation is to visit more and more neighbors. Moreover, in a hyperbolic framework, the aim is to rebuild a non oscillatory interpolation [START_REF] Harten | Some results on uniformly high-order accurate essentially nonoscillatory schemes[END_REF], which complicates the problem. This method is not compact, and therefore is not well suited to parallel computing.

• Continuous Galerkin method

In this method, the solution is approximated by piecewise continuous polynomial functions. The numerical scheme is then obtained by writing the L 2 orthogonality between the approximation basis and the equation projected on the approximation space. It is well adapted to some problems such as purely parabolic or elliptic problems. For Euler equations, and more generally for hyperbolic systems (even linear), this method is known to be unstable. This method can be stabilized, for example with the SUPG method [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]. Nevertheless, this method depends on parameters that can be hard to tune.

• Residual Distribution schemes

The development of fluctuation schemes began with the early work of Phil Roe [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF].

Residual distribution is a weighted residual approach in which local discrete equations are derived as a sum of elemental contributions proportional to element integrals of the equations (the cell residuals) via matrix weights. The basics of the method are thoroughly discussed in [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. As in continuous and discontinuous finite element methods, the key toward high accuracy is the use of a high order polynomial representation of the unknowns in the computation of the cell residuals [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF]. Although this approach has shown great potential in steady applications [START_REF] Abgrall | Contruction of very high order residual distribution for steady inviscid flow problems on hybrid unstructured mesjes[END_REF][START_REF] Abgrall | Some examples of high order simulations of inviscid flows on unstructured hybrid meshes by residual distribution schemes[END_REF], its current state of the art [START_REF] Abgrall | A review of residual distribution schemes for hyperbolic and parabolic problems : the july 2010 state of the art[END_REF] shows that further work is needed to bring the method to the level of maturity of more popular techniques, such as Discontinuous Galerkin.

• Discontinuous Galerkin method The development of discontinuous Galerkin for nonlinear hyperbolic equations began in [START_REF] Chavent | The local projection P 0 P 1 -discontinuous-Galerkin finite element method for scalar conservation laws[END_REF]. Its stabilization for flows with shocks was developed in the 90's, mainly by Cockburn and Shu (see [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] for a review). In the same time, a solution for dealing with Navier-Stokes equations was proposed in [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF]. In spite of its cost (it has much more degrees of freedom than classical continuous finite elements methods), it is attractive because it is naturally L 2 stable for linear problems, because a cell entropy inequality can be proven [START_REF] Shan | On a cell entropy inequality for discontinuous Galerkin methods[END_REF], and also because it has a compact stencil so that it has a good behavior in parallel environment [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF]. The success of these methods lies in their flexibility thanks to their high degree of locality. These properties make the DG method well suited to parallel computing, hp-refinement, hp-multigrid, unstructured meshes, the weak application of boundary conditions, etc.

Inria

The development of different high order methods for aerospace application was the topic of the European project ADIGMA, and we refer to [START_REF] Kroll | ADIGMA -A European Initiative on the Development of Adaptive Higher-Order Variational Methods for Aerospace Applications[END_REF] for recent developments on this topic.

1 Discontinuous Galerkin methods and libraries

The Euler model

Let Ω ⊂ R d be a bounded domain where d is the space dimension and consider the following problem

∂ t u + ∇ • f (u) = 0, in Ω × (0, ∞) , (1)
with initial condition u(•, 0) = u 0 (•) in Ω and appropriate boundary conditions prescribed on ∂Ω. The vector u = (ρ, ρv, ρE) represents the conservative variables with ρ the density, v ∈ R d the velocity vector and E = ε(p, ρ) + v 2 /2 where ε is the specific internal energy. Here, we suppose that the fluid follows the perfect gas equation

ε(p, ρ) = p (γ -1)ρ ,
where p denotes the pressure and γ is the ratio of specific heats. The nonlinear convective fluxes in (1) are defined by

f (u) =      ρv ρvv + pI (ρE + p)v      . ( 2
)
The problem (1) is hyperbolic provided the conservative variable vector takes values in the set of admissible states

Ψ = {u ∈ R d+2 : ρ > 0, E - v 2 2 > 0} . (3)

Runge-Kutta discontinuous Galerkin formulation

The discontinuous Galerkin method is a finite element method in which the weak formulation is projected on a space of discontinuous piecewise polynomials of the problem [START_REF] Abgrall | A review of residual distribution schemes for hyperbolic and parabolic problems : the july 2010 state of the art[END_REF]. The domain Ω is partitioned into a shape-regular mesh Ω h consisting of nonoverlapping and nonempty elements κ of characteristic size h := min{h κ , κ ∈ Ω h } where h κ is a d-dimensional measure of κ. We define the sets E i and E b of interior and boundary faces in Ω h , respectively. We look for approximate solutions in the function space of discontinuous polynomials

V p h = {φ ∈ L 2 (Ω h ) : φ| κ • F -1 κ ∈ P p (κ)
, ∀κ ∈ Ω h }, where P p (κ) denotes the polynomial space associated to the reference element κ corresponding to the element κ. Each physical element κ is the image of one of the following reference shapes κ through the mapping F κ : simplex (line, triangle or tetrahedron), tensor elements (quadrangle, hexaedron), prism or pyramid. The space P p (κ) might be composed of the set of polynomials with degree lower or equal to p in the case of simplex, or of tensor product of one dimensional polynomials with degree lower or equal to p in the case of tensor elements, or of a tensor of one dimensional basis and two dimensional basis on a triangle if κ is a prism, or a conical product of one dimensional basis and two dimensional basis on a quadrangle if κ is a pyramid. The numerical solution of equation ( 1) is sought under the form

v h + v h - n e κ + κ -
u h (x, t) = Nκ l=1 φ l κ (x)U l κ (t), ∀x ∈ κ, κ ∈ Ω h , ∀t ≥ 0 , (4)
where (U l κ ) 1≤l≤Nκ are the degrees of freedom in the element κ. The semi-discrete form of equation ( 1) reads: find u h in [V p h ] d+2 such that for all v h in V p h we have

Ω h v h ∂ t u h dx + B h (u h , v h ) = 0 . (5)
Hereafter, we will use the notation [[φ]] = φ + -φ -which denotes the jump operator defined for a given interface e ∈ E i . Here, φ + and φ -are the traces of any quantity φ on the interface e taken from within the interior of the element κ + and the interior of the neighboring element κ -, respectively (see Figure 1).

The space discretization in equation ( 5) can then be written as

M κ u n+1 h -u n h δt + B h (u n h , v h ) = 0 . (6)
for the explicit Euler time stepping, the extension to other explicit time stepping being straightforward. In (6) M κ denotes the local mass matrix defined as

∀(i, j) ∈ [0; N κ ] × [0; N κ ] M (i,j) κ = κ φ i κ (x)φ j κ (x) dx ,
and the space discretization operator B h is defined by

B h (u h , v h ) = - Ω h f (u h ) • ∇v h dV + Ei [[v h ]] f (u + h , u - h , n)dS + E b v + h f (u b (u + h , n)) • ndS , (7)

Inria

where n denotes the unit outward normal vector to an element κ + (see Figure 1) and u b is an appropriate operator which allows to impose boundary conditions on E b . The numerical flux f may be chosen to be any monotonic Lipschitz function satisfying consistency, conservativity and entropy dissipativity properties (see [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] for instance). Throughout this study, the semi-discrete equation ( 5) is advanced in time by means of an explicit treatment with third-order strong stability preserving Runge-Kutta methods [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF][START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF].

Algorithm

As an explicit Runge-Kutta time stepping was chosen, the algorithm consists in computing the spatial residual, and then in updating the intermediate (or the final) steps of the Runge-Kutta method by inverting the mass matrix. In the discontinuous Galerkin method, the mass matrix is block diagonal. It can actually be diagonal provided an orthogonal basis is used on each element. We point out that this step of the method does not require any communication in a distributed memory environment if all elements are on the same process (which is always the case). We now dwell on the different loops that must be made for computing the local spatial residual defined on [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF]. We detail the computations needed for the elements loop, the other loops being similar.

We are interested in computing

κ f (u h ) • ∇v h dV
for all basis function v h of κ. Using the definition of the basis given in Section 1.

2 κ f (u h ) • ∇φ i κ dV = κ f ( Nκ l=1 φ l κ (x)U l κ ) • ∇φ i κ dV .
As in Section 1.2, we denote by F κ the map from κ to κ. In the integral, we change the variable x by x = F κ (x), so that

κ f Nκ l=1 φ l κ (x)U l κ • ∇φ i κ dV = κf Nκ l=1 φ l κ • F κ (x)U l κ • DF -1 κ ∇ φi κ |det DF κ | dx = κf Nκ l=1 φl κ(x)U l κ • DF -1 κ ∇ φi κ |det DF κ | dx .
An approximated quadrature formula is used for computing this integral. We denote by xα the quadrature points and by ω α the weights.

κ f Nκ l=1 φ l κ (x)U l κ •∇φ i κ dV ≈ α ω α f Nκ l=1 φl κ(x α )U l κ •DF -1 κ (x α )∇ φi κ (x α ) |det DF κ (x α )|

Remarks on quadrature formulas

For linear problems, quadrature formulas can be chosen for being exact. For nonlinear problems, [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] suggests, for an approximation of degree p, to take a (2p)th order formula for cells, and a (2p + 1)th order for faces. For hypercube shapes, the optimal set of points (i.e. the one ensuring the highest degree with a given number of points) is the Gauss set of points. For other shapes, the optimal set of points is often unknown. A systematic way of deriving quadrature formula on simplexes consists in using the image of Gauss points by Dubiner's map [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF]. As far as we know, this is what is done in the quadrature formula proposed in [START_REF] Šolín | Higher-order finite element methods[END_REF]. However, the quadrature formulas obtained are not symmetric, and the number of points might be greater than other sets obtained by optimization procedures, see e.g. [START_REF] Zhang | A set of symmetric quadrature rules on triangles and tetrahedra[END_REF], and [START_REF] Stroud | Approximate calculation of multiple integrals[END_REF][START_REF] Cools | An encyclopaedia of cubature formulas[END_REF] for a list of quadrature formulas on simplexes.

The AGHORA and AEROSOL libraries

The AGHORA and AEROSOL libraries are two high order finite elements libraries that are respectively developed within ONERA and INRIA Bordeaux Sud Ouest (BACCHUS and CAGIRE teams). In this section, we give some details on these libraries.

The AGHORA library

The development of the AGHORA library began in 2012 with the PRF 1 of the same name. It is a high order (arbitrary order) finite elements library based on discontinuous elements. The models solved are the 3d compressible Euler or Navier-Stokes and RANS equations (with BR2 formulation for diffusive fluxes, see [START_REF] Bassi | Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k -ω turbulence model equations[END_REF]) based on straightsided and curved tetrahedra, hexahedra and prisms. It is developed in FORTRAN 95, and its development is based on previous experiences on two dimensional codes, that were developed for example within the ADIGMA project.

The AEROSOL library

The development of the AEROSOL library began at the end of 2010 with the PhD of Damien Genet. It became a shared project between the teams BACCHUS and CAGIRE in 2011. It is a library that aims at using tools that are developed mainly within Inria teams working on high performance computing at the Bordeaux center, see Figure 2.

Under the bonnet

AEROSOL is a high order finite elements library based on both continuous and discontinuous elements on hybrid meshes involving triangles and quadrangles in two dimensions and tetrahedra, hexaedra and prisms in three dimensions. More precisely, the finite element classes are generated until the fourth order polynomial approximation. Currently, it is possible to solve simple problems with the continuous Galerkin method (Laplace equation, SUPG stabilized advection equations), with the discontinous Galerkin method (Hyperbolic systems of first order) and with residual distribution schemes (scalar hyperbolic equations). It is written in C++, and its development started nearly from scratch as far as the general structure of the code is concerned. It strongly depends on the PAMPA library (see next section for details). It is linked with external linear solvers (up to now, PETSC 2 and MUMPS 3 ). It is about to use the STARPU 4 task scheduler also being developed at INRIA. The choice of C++ allows for a good flexibility in terms of models and equations of state. Currently, the following models It is about to be linked with the STARPU library, which is a task scheduling library for hybrid architectures. can be used: scalar advection, waves in a first order formulation, nonlinear scalar hyperbolic equation and Euler model with an abstract equation of state (currently: perfect gas and stiffened gas equation of state). It is working on linear elements, but the level of abstraction is sufficient for taking into account curved elements.

The PAMPA library

The PAMPA middleware library aims at abstracting mesh handling operations on distributed memory environments. It relieves solver writers from the tedious and error prone task of writing service routines for mesh handling, data communication and exchange, remeshing, and data redistribution. It is based on a distributed graph data structure that represents meshes as a set of entities (elements, faces, edges, nodes, etc.), linked by relations (that is, computation dependencies).

Given a numerical method based on a mesh, the user shall define an entity graph containing all the entities holding an unknown. For example, in a discontinuous Galerkin formulation, all the unknowns are located on the cells of the mesh, whereas in high order continuous finite elements, the unknowns might lie not only on elements, but also on points, edges and faces. Given the entity graph, PAMPA is able to compute a balanced partition of the unknowns (see Figure 4) and to compute adequate overlaps for data exchange across processors (see Figure 5). PAMPA handles mesh memory allocation, so that it can release memory after remeshing and/or mesh redistribution, as well as communication (roughly speaking, the user has almost no MPI to write in his code). PAMPA provides methods for iterating on mesh entities (e.g. on all local or boundary elements, on faces of an element, etc.), which eases the writing of numerical schemes based on finite elements methods. Last, PAMPA is also able to perform mesh adaptation in parallel, provided a sequential mesh adaptation software is linked to it; see Figure 3 for more details. Overview of the use of PAMPA and of its interaction with other software. Solver writers (top box) rely on PAMPA for all mesh structure and data handling. PAMPA strongly depends on PT-SCOTCH for partitioning and redistribution. In order to perform (optional) dynamic mesh refinement with PAMPA, solver writers have to provide a sequential mesh refinement module (bottom box) that handles their types of elements, as well as a sequential quality measurement metric (rightmost box) to tell PAMPA where to perform remeshing. Remeshing is performed in parallel, after which the refined mesh is redistributed by PAMPA so as to re-balance computation load.

Comparison of implementations 2.3.1 How to take into account the geometry

The main difference concerning implementations regards the strategy for taking into account geometries of cells in the different loops. For example, if we are concerned with the element loop, we have to compute

κ f Nκ l=1 φ l κ (x)U l κ •∇φ i κ dV ≈ α ω α f Nκ l=1 φl κ(x α )U l κ •DF -1 κ ∇ φi κ (x α ) |det DF κ | .
In this formula, some of the terms do not depend on the geometry of the cell: ω α , φl κ(x α ), and ∇ φi κ (x α ), whereas the following terms depend on the geometry: DF -1 κ and |det DF κ |. The strategy in AGHORA was to store these geometrical dependent terms. This is very memory costly, because one needs to store it on all quadrature points (if elements are not linear, the geometrical terms are not constant in one given element), but it is often considered as paying off, as their evaluation is also costly. In the AEROSOL library, the only item stored is a pointer to the function that computes these geometrical terms. It is called for all cells at each computation step.

Computation of basis

In the same manner, the way to compute finite element basis is more generic in the AGHORA library. In the AEROSOL library, everything is based on the reference element. An orthogonal basis on the reference element is such that

κ φi κ(x) φj κ(x) dx = δ i,j κ φi κ(x) 2 dx ,
where δ i,j is the Kronecker symbol. If we consider an element κ with associated reference elements κ, i.e. such that F κ (κ) = κ, and if we consider the basis φi

κ • F -1 κ , then κ φ i κ (x)φ j κ (x) dx = κφ i κ • F κ (x)φ j κ • F κ (x) |det F κ | dx = κ φj κ(x) |det F κ | dx . If F κ is linear, then |det F κ | is constant, and the basis φi κ • F -1 κ is orthogonal on κ. But if F κ is nonlinear, for example if κ is a hexaedra which is not a parallelepiped, or if κ is a curved simplex, then the basis φi κ • F -1
κ may not be orthogonal. The library AGHORA is designed such as the finite element basis is always orthogonal.

For the sake of clarity, we introduce the construction of the orthonormal basis in the 1D case. The generalization to multi-space dimensions is straightforward. We start with a Taylor expansion of the numerical solution about x κ the barycentre of the element κ:

u h (x, t) = Nκ l=1 φl κ (x) Ûl κ (t), ∀ x ∈ κ, ∀t ≥ 0 , (8)
where

Inria φ1 κ (x) = 1 , φl κ (x) = (x -x κ ) l-1 -(x -x κ ) l-1 κ (l -1)! , ∀2 ≤ l ≤ N κ , x ∈ κ . (9)
Functions φl κ , for 2 ≤ l ≤ N κ are centered monomials with zero mean value and • κ denotes the average operator over the element. The degrees of freedom in [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] are defined by

Û1 κ (t) = u h κ , Ûl κ (t) = ∂ l-1 u h ∂x l-1 x=xκ , ∀2 ≤ l ≤ N κ , t ≥ 0 . (10)
Functions ( φl κ ) 1≤l≤Nκ are not orthogonal with respect to the inner product defined over the element κ and lead to a non-diagonal and ill-conditioned mass matrix M κ even for p ≥ 1 if d ≥ 2. Their inversion for the time integration therefore requires extra computational costs and convergence properties of the method deteriorate. For general meshes, it is however possible to construct an orthonormal basis (φ l κ ) 1≤l≤Nκ in an arbitrary element κ by applying a Gram-Schmidt orthonormalization to the initial basis ( φl κ ) 1≤l≤Nκ . The new basis satisfies the following orthonormality property

(φ k κ , φ l κ ) κ = κ φ k κ (x)φ l κ (x)dV = |κ|δ kl , ∀1 ≤ k, l ≤ N κ , (11)
and the mass matrix reduces to a diagonal matrix of the form M κ = |κ|I. We refer to [START_REF] Bassi | On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations[END_REF][START_REF] Remacle | An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems[END_REF] for details on this procedure. As a consequence of the orthonormalization algorithm, both basis in expansions ( 4) and ( 8) are related by the following relation

φl κ (x) = l k=1 r lk κ φ k κ (x) , ∀1 ≤ l ≤ N κ , ∀x ∈ κ , (12)
where r lk κ = ( φl κ , φ k κ ) κ , for 1 ≤ k ≤ l -1, denotes the inner product [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF], and

(r ll κ ) 2 = ( φl κ , φl κ ) κ - l-1 k=1 ( φl κ , φ k κ ) 2 κ .

Communication handling

AGHORA performs data exchange by means of point-to-point communications, whereas AEROSOL uses collective communications using the overlap data exchange routines provided by PAMPA. If POSIX THREADs are available, PAMPA can perform such collective communications asynchronously on a specific thread, thus overlapping communication with computation. We were not yet able to fully use this feature, as most highspeed MPI implementations do not support the MPI_THREAD_MULTIPLE model.

Comparisons

Numerical tests

Due to the fact that both the AEROSOL and AGHORA libraries are not yet mature, a simple test was chosen for comparing libraries. In AGHORA, only three dimensional computations are possible. One and two dimensional computations are also possible, but by extruding one layer of cells in the z direction (and also in the y direction in one dimension). In AEROSOL, true one and two dimensional computations can be considered, but geometrical functions, finite elements functions, and mesh reading were not yet available for three dimensional shapes. For performing fair comparisons, the AGHORA point of view was adopted, and one layer hexaedral meshes were used for doing a two dimensional test. Thus, AEROSOL was extended to three dimensions.

Yee vortex

We consider the convection of an isentropic vortex in a 2D uniform and inviscid flow [START_REF] Yee | Low-dissipative high-order shock-capturing methods using characteristic-based filters[END_REF] with conditions ρ ∞ = 1, u ∞ = e x and T ∞ = 1. The domain is the unit square Ω = [0, 1] 2 with periodic boundary conditions. The initial condition consists in a perturbation of the uniform flow which reads in primitive variables

ρ(x, 0) = 1 -(γ -1) M∞Mvrc 2 exp(1 -r 2 r 2 c ) 2 1 γ-1 , u(x, 0) = 1 -M v y exp(1 -r 2 r 2 c ) , v(x, 0) = M v x exp(1 -r 2 r 2 c ) , p(x, 0) = 1 γM 2 ∞ ρ(x, 0) γ ,
where r 2 = (x -x c ) 2 + (y -y c ) 2 denotes the distance to the vortex centre (x c , y c ) , r c and M v are the radius and strength of the vortex, and M ∞ is the Mach number of the freestream flow. The exact solution of this problem is a pure convection of the vortex at velocity u ∞ .

Numerical results are obtained for physical parameters M ∞ = 0.5, M v = 4, r c = 0.1, and a final time T = 1.

As the periodic boundary conditions were not yet available in the AEROSOL library, the test was slightly modified as follows: instead of adding a mean flow to the vortex, the initial state is taken as the state at infinity, without any velocity. On the boundaries, the values of the vortex are weakly applied. From a computational point of view, the work load is nearly the same, except for the communications that are needed for periodic boundary conditions in the case of the unsteady vortex.

Table 1 presents a weak scalability analysis where we evaluate the elapsed time observed for the global computation and for a fixed number of 24 × 24 hexaedral elements per computing core. Results are shown as a function of the number of cores. The relative increase, E, of the elapsed time with respect to the single core is also indicated shown, on Figure 6. Note that the memory requirement increases with the polynomial degree: the number of degrees of freedom per core is 4608, 15552, and 36864 for p = 1, 2, and 3, respectively. We observe that the solver scales nearly perfectly for p = 4, while results deteriorate for lower polynomial degrees p = 1 and 2. The relative part of communications compared to the numerical scheme decreases as the polynomial degree increases : from 1.1% for p =1 to 0.1% for p=4, as we can see in Table 2. As a consequence, the high frequency of communications between processes for a large number of cores becomes insensitive for large p-values. AEROSOL has a stronger work load for boundary sides, because in AGHORA, the boundary sides of the top and the bottom are ignored in the boundary side loop, whereas in AEROSOL, a freestream boundary condition was used. All the results show good scalability behavior. The least efficient is AGHORA for p = 1, and the most efficient is AGHORA for p = 4. Nevertheless, if we consider the results given in Table 1, we see that AEROSOL is much more costly if execution time is considered. This explains that the weak scalability is better, as AEROSOL has more computations for hiding communications.

Inria

We observe in Table 1 that AGHORA has lower execution time for a same number of RK substeps. Using an analyzer, we found the following reasons for the lower efficiency of AEROSOL • First, AEROSOL is a finite element code that can both use continuous and discontinuous finite elements. In continuous finite elements, a degree of freedom is shared by many elements, if it is located on edges, faces, or points. Consequently, a connectivity table must be built for knowing which degrees of freedom is shared by which element. In AEROSOL, this connectivity table is never stocked, and its computation is general for continuous and discontinuous finite elements (thus it was not optimized for discontinuous finite elements). We found that the computation of this connectivity table is nearly three times more costly that the computation of the local residual. We are currently developing a version in which the connectivity table is stocked.

• Second, in AEROSOL, the geometry (e.g. local Jacobian of elements) is never stocked. This means that at each time step, and on all quadrature points, all the needed geometry is computed. In the computation of the local residual, we evaluated that 35% of the time elapsed in the computation of the local residual is for computing geometrical functions. We will develop in the near future a version in which the geometry is stocked.

• Last, we did not use any optimization option for compiling AEROSOL.

Table 2 gives a detailed analysis of the relative costs of different stages of the numerical algorithm. Results are given for one core. The implementation of surface and volume integrals for fluxes clearly represents the most expensive stage and its relative cost increases with p. This result is in agreement with the high scalability observed for high p-values as local operations strongly dominate communications.

Last, we want to compare the results obtained with different MPI libraries; however, AEROSOL needs an MPI library for which the MPI_THREAD_MULTIPLE is supported (i.e. if the process is multithreaded, multiple threads may call MPI at once with no restrictions), and such a functionality is only available with MVAPICH2 on AVAKAS. That is why we show the results of this comparison only for AGHORA and not for AEROSOL, see Table 3. Note that an important part of this execution time can be due to the initilization of the MPI universe, see Table 4.

Conclusion

This project was the opportunity for us to compare the efficiency of the libraries AEROSOL and AGHORA. A weak scalability analysis was performed for comparing the two codes. We found that the building of the local connectivity and the geometry was very costly, and can explain a gap between the performances of the two codes. Concerning the weak scalability, we found that both of the codes have a good behavior, at least until 2048 cores.

Further than the results presented here, this project allowed us to have discussions on the way to develop and to factorize the code, that are difficult to account for here.

Figure 1 :

 1 Figure 1: Inner and exterior elements κ + and κ -and definition of traces v ± h on the interface e and of the unit outward normal vector n.

Figure 2 :

 2 Figure 2: Structure of the AEROSOL code. It uses the PAMPA library for memory handling, for mesh partitioning, and for abstracting the MPI layer. It uses external linear solver (currently: PETSC and MUMPS).It is about to be linked with the STARPU library, which is a task scheduling library for hybrid architectures.

  Figure3: Overview of the use of PAMPA and of its interaction with other software. Solver writers (top box) rely on PAMPA for all mesh structure and data handling. PAMPA strongly depends on PT-SCOTCH for partitioning and redistribution. In order to perform (optional) dynamic mesh refinement with PAMPA, solver writers have to provide a sequential mesh refinement module (bottom box) that handles their types of elements, as well as a sequential quality measurement metric (rightmost box) to tell PAMPA where to perform remeshing. Remeshing is performed in parallel, after which the refined mesh is redistributed by PAMPA so as to re-balance computation load.

Figure 4 :

 4 Figure 4: Example of mesh redistribution performed by PAMPA. Initially, each process reads a piece of a sequential mesh file (left picture), which may lead to poor data locality. Once the entity graph is given to PAMPA, it is able to compute a mesh partition that balances computation across processes and minimizes communication load (right picture).

Figure 5 :

 5 Figure 5: Examples of overlaps that can be computed by PAMPA according to the requirements of the numerical schemes. Fig. (a) represents the redistribution of the cells of a mesh across four processors, as yielded by PAMPA. Fig. (b) shows the overlap corresponding to a P 1 continuous finite element method (node-based neighbors), Fig. (c) shows the overlap for a discontinuous Galerkin or cell-centered finite volume method (face-based neighbors), and Fig. (d) shows the overlap for a high order finite volume scheme. Every processor of a given color stores its local data (Fig. (a)) plus the overlap cells of the same color.

Figure 6 :

 6 Figure 6: Comparison of the weak scalability obtained for AEROSOL and AGHORA with the library MVA-PICH2.All the results show good scalability behavior. The least efficient is AGHORA for p = 1, and the most efficient is AGHORA for p = 4. Nevertheless, if we consider the results given in Table1, we see that AEROSOL is much more costly if execution time is considered. This explains that the weak scalability is better, as AEROSOL has more computations for hiding communications.

1

  Projet de Recherche Fédérateur 2 www.mcs.anl.gov/petsc/ 3 graal.ens-lyon.fr/MUMPS/ 4 starpu.gforge.inria.fr/

	Aerosol
	PAMPA Linear Solver StarPU
	PT-Scotch
	Inria

Table 1 :

 1 Yee vortex problem: time/proc. [s] with Aerosol and Aghora for 3000 Runge-Kutta substeps. Aerosol 202.52 205.32 209.91 216.68 222.69 224.01 225.52 Aghora 94.12 96.90 100.08 100.72 103.28 112.91 116.75 Aerosol 5938.6 5923.2 6047.1 6307.7 6504.6 6501.0 6640.5 p = 4 Aghora 1401.5 1421.3 1430.2 1435.5 1439.1 1451.4 1370.8

	# cores	1	4	16	64	256	1024	2048
	p = 0 p = 1 Aghora	12.67	14.02	15.69	15.88	16.42	18.43	17.97
	Aerosol 977.3	977.1	999.3 1035.6 1074.5 1074.3 1083.9
	p = 2

Table 2 :

 2 Yee vortex problem: relative costs in percent of different stages of the numerical algorithm per physical time step.

	p	0	1		2		4
		Aerosol Aerosol Aghora Aerosol Aghora Aghora
	boundary surface integral	31.1	25.0	4.2	16.5	3.2	2.2
	internal surface integral	51.3	36.4	57.5	21.2	46.8	34.6
	volume integral	17.3	37.6	36.5	61.4	49.3	63.0
	mass matrix inversion	0.3	1.0	1.8	0.8	0.7	0.2
							Inria

Table 3 :

 3 Comparison of execution time by core for a polynomial degree equal to 2 (9000 RK substeps). We observe that the most efficient execution is obtained for the OpenMPI library on AVAKAS.

	# core	256 1024 2048
	IntelMPI	15	122 5001
	OpenMPI	05	010 0020
	MVAPICH2 33	122 0536

Table 4 :

 4 Average elapsed time (s) required for initialisation of MPI universe.

RR n°8200

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements Part of the computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain, on the cluster AVAKAS) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. Part of the computations have been performed at the Mésocentre d'Aix-Marseille Université.


Ismail Assayad
email: ismail.assayad@ensem.ac.maalain.girault@inrialpes.frhamoudi.kalla@univ-batna.dz

Alain Girault

Hamoudi Kalla

ENSEM Rtse Team

Scheduling of Real-Time Embedded Systems under Reliability and Power Constraints

Keywords: Distributed real-time systems, safety-critical systems, reliability, multi-criteria scheduling, heterogeneous systems, active software replication

In this paper, we present a new tri-criteria scheduling heuristic for scheduling data-flow graphs of operations onto parallel heterogeneous architectures according to three criteria: first the minimization of the schedule length crucial for real-time systems, second the maximization of the system reliability crucial for dependable systems, and third minimizing energy consumption crucial for autonomous systems. The proposed algorithm is a list scheduling heuristics, It uses the active replication of operations to improve the reliability and the dynamic voltage scaling to minimize the energy consumption.

I. Introduction

Distributed systems are being increasingly used in critical real-time applications, such as avionics, air traffic control, autopilot systems, and nuclear plant control, in which the consequences of missing a tasks deadline may cause catastrophic loss of money, time, or even human life. This is why such systems require a high reliability.

Here, reliability is defined as the probability that none of the system components will fail while processing.

We present a scheduling heuristics that, from a given software application graph and a given multiprocessor architecture, produces a static multiprocessor schedule that optimizes three criteria: its length (crucial for real-time systems), its reliability (crucial for dependable systems), and its power consumption (crucial for autonomous systems). Let us address the issues raised by multicriteria optimization. Figure 1(a) illustrates the particular case of two criteria to be minimized. Each point x 1 to x 7 represents a solution, that is, a different tradeoff between the Z1 and Z2 criteria: the points x 1 , x 2 , x 3 , x 4 , and x 5 are Pareto optima [13]; the points x 1 and x 5 are weak Ismail Assayad is with ENSEM (RTSE team). Partially funded by DEV-RS 2011/14 CNRST project.

Second criterion Z 2 x 3 x 2 x 6 x 1 x 7 x 4 x 5 Z 1 (x 6 ) First criterion Z 1 Z 2 (x 6 )
Alain Girault is with INRIA and Grenoble University (POP ART projectteam and LIG laboratory), France.

Hamoudi Kalla is with the University of Batna (REDS team), Algeria.

optima while the points x 2 , x 3 , and x 4 are strong optima. The set of all Pareto optima is called the Pareto curve.

It is fundamental to understand that no single solution among the points x 2 , x 3 , and x 4 (the strong Pareto optima) can be said, a priori, to be the best one. Indeed, those three solutions are non-comparable, so choosing among them can only be done by the user, depending on the precise requirements of his/her application.

The main contribution of this paper is TSH, the first tricriteria scheduling heuristics able to produce, starting from an application algorithm graph and an architecture model, a set of Pareto solutions in the space (length, reliability, power), and taking into account the relationship between the reliability and the voltage level. In addition, the power consumption of both operations and data-dependencies are considered in this paper and the soundness of TSH is demonstrated.

II. Related work

Many solutions exist in the literature to optimize the schedule length and the energy consumption (e.g., [10]), or to optimize the schedule length and the reliability (e.g., [5]), but very few tackle the problem of optimizing the three criteria (length,reliability,energy). The closest to our work are [11], [14].

Zhu et al. have studied the impact of the supply voltage on the failure rate [14], in a passive redundancy framework (primary backup approach). They use DFS and DVS to lower the energy consumption and they study the tradeoff between the energy consumption and the performability (defined as the probability of finishing the application correctly within its deadline in the presence of faults).However, their input problem is not a multiprocessor scheduling one since they study the system as a single monolithic operation executed on a single processor.

Pop et al. have addressed the tricriteria optimization problem (length,reliability,energy) on an heterogeneous architecture with a reliable bus [11]. Both the length and the reliability are taken as a constraint, respectively with a given upper and lower bound. However, it is assumed that the user will specify the number of processor failures to be tolerated (with re-execution or passive replication) in order to satisfy the desired reliability constraint.

III. Models

A. Application algorithm graph

Most embedded real-time systems are reactive, and therefore consist of some algorithm executed periodically, triggered by a periodic execution clock. Our model is therefore that of an application algorithm graph Alg which is repeated infinitely. Alg is an acyclic oriented graph (X , D). Its nodes (the set X ) are software blocks called operations. Each arc of Alg (the set D) is a data-dependency between two operations. If X⊲Y is a data-dependency, then X is a predecessor operation of operation Y , while Y is a successor operation of X. Operation X is also called the source of the data-dependency X ⊲ Y ,

Fig. 1 .

 1 Fig. 1. Pareto optima and Pareto curve for a bicriteria minimization problem.