Jean Marie Couteyen Carpaye

Jean Roman

Pierre Brenner

Towards an efficient Task-based Parallelization over a Runtime System of an Explicit Finite-Volume CFD Code with Adaptive Time Stepping

FLUSEPA 1 is an advanced simulation tool which performs a large panel of aerodynamic studies. It is the unstructured finite-volume solver developed by Airbus Defence & Space 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 a first study of a task-based parallelization of the solver part of FLUSEPA using the runtime system StarPU and combining up to three levels of parallelism. We validate our solution on the simulation of a take-off blast wave propagation for Ariane 5 launcher.

I. 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 1 Registered trademark in France No. 134009261 different kinds of architectures (performance portability issue [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF], [START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF]).

The aerospace industry faces a lot of complex problems for which actual experiences are not doable: e.g. take-off blast wave propagation, 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[END_REF] of the aerodynamic solver of the FLUSEPA code [START_REF] Brenner | Simulation d'interactions aérodynamiques instationnaires autour de plusieurs corps en mouvement relatif[END_REF] which is an MPMD MPI/OpenMP code. This code is able to compute stage separation. Each body is meshed separately and can move using an ALE formulation [START_REF] Hirt | An arbitrary lagrangianeulerian computing method for all flow speeds[END_REF], leading to mesh intersection computations. The code is also well adapted for unsteady computations. Figure 1 shows a separation computation: each booster and the main stage are meshed separately. The aerodynamic solver of FLUSEPA 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] Lhner | 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).

The paper is organized as follows. Section II briefly presents the main computational methods used in FLUSEPA and the existing MPI/OpenMP code. Section III presents the runtime system StarPU and Section IV focuses on the task parallelization of the aerodynamic solver. Section V presents an experimental study from an industrial test case, the takeoff blast wave propagation. Finally, Section VI gives some conclusions and perspectives for this work.

II. MAIN NUMERICAL METHODS USED IN FLUSEPA

The two main computation steps performed in FLUSEPA, 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. This general computational framework is summed up in Algorithm 1.

Algorithm 1 General iteration

Intersection computation 6: end if

A. Finite-volume aerodynamic solver

Finite-volume method is a discretization technique where the integral formulations of the conservation laws are discretized directly in the physical space. The finite-volume method is naturally conservative. We look for the numerical solution of the compressible Navier-Stokes equations in Reynolds-Averaged form 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. See [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF] for more details.

B. 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 Predictor: end for 36: end 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 θ = 2, τ will take successively the following values : 2,0,1,0. Figure 2 shows how level τ evolves after each subiteration in this case; the evolving levels τ are in red. The physical time reached after one iteration with θ = 2 is equivalent to 4 iterations with a global time step.

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.

C. 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].

We focus here on the aerodynamic solver computations. 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 3 illustrates this spatial decomposition with 2 subdomains: for the red subdomain, the dark red part corresponds to the border cells, the light red part to inner cells, the dark green cells are the ghost cells of this red subdomain, and finally bold black 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 computationcommunication 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. This 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.

III. 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[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[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 pop 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.

A. 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 workercontexts and insert tasks in this main-context.

B. Distributed parallelism with StarPU

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.

It is also possible to describe communications explicitly using starpu_mpi_isend_detached / starpu_mpi_irecv_detached. When using those primitives, communications are consistent with the computation dependencies inferred by task insertion inside the nodes.

IV. MPI + TASK DESIGN OF THE AERODYNAMIC SOLVER A. 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 SCOTCH 2 [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 II-A). The different components of a CE are shown at Figure 3: there are two CEs, a red one and a green one. Border cells (dark red and dark green) and inner cells (light red and light green) of the CEs are distinguished. In a CE, faces are classified into different categories: faces between inner cells (black), faces between border cells (white), faces between border/inner cells (bold light blue). Between CEs, we also have inter-CE faces (bold black). Because we have a unstructured mesh, all these topological informations are precomputed for each CE after the partitioning step.

B. Task generation

In order to generate tasks, we mainly use the different data types we described previously and exploit the associated 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 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.

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 I which gives the proportion of cells in each temporal level and the proportion of associated computation for the test case we use in Section V. The large majority of cells are of the higher temporal level (τ 4 ). 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. 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.

C. 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.

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 Fig. 4: DAG for a computation with 6 CEs and θ = 2.

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.

V. 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 5). 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 6: the comeback of the wave can be seen from the 4th picture.

We consider a case with θ = 4. The cell distribution in temporal levels and the theoretical computation load associated to each level is shown in Table I.

We focus on one iteration of the aerodynamic solver and we study the impact of several parameters: the priority strategy for the scheduling of tasks, the number of CEs, the number of parallel workers. We start by a study in shared memory and then we will evaluate the distributed version of our new solution.

A. Shared memory study 1) 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 7. 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 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 green (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 priority they have in the computation with priorities. There are four priority levels: from dark blue (high priority task) to light blue (low priority task). The last operation (line 15 of Algorithm 3) is in light green.

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 low priority tasks.

The evolution of ready tasks in the solver is shown at Figure 9. 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 shared memory study, 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. 2) 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 worker)×(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 10, 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.

B. A distributed memory version

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 11 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 II. 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 communicationcomputation overlapping.

Figure 11d shows the elapsed time for the previous 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.

VI. CONCLUSION AND PERSPECTIVES

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

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. However, for our current fixed-sized problem, the overhead becomes too high to obtain a gain when we use too much nodes. As we expect to grow the size of our industrial test-case, this is not dramatic. So this result is very promising and the gain of our distributed task-based solution should be more important with a larger size test case.

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.

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.

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].

Concerning the distributed version of the code, we have to test the implementation with a bigger test case. If we manage to remove the hard synchronization we have at the end of iterations, we 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.

Fig. 1 :

 1 Fig. 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 θ + 1 do for tmp = 1 to θ + 1 do 7:if (mod(subiteration -1, 2 tmp ) == 0) then

Fig. 2 :

 2 Fig. 2: Time reached after each subiteration for each level τ (θ = 2; 4 subiterations).

Fig. 3 :

 3 Fig. 3: Illustration of two CEs (a red one and a green one).

Algorithm 3 5 :

 35 2 https://gforge.inria.fr/projects/scotch/ Task insertion for the Aerodynamic Solver 1: Insert tasks for time step computation 2: Insert tasks for classification of cells in temporal levels 3: Wait all tasks 4: Compute hints (scheduling and packing algorithm) Temporal adaptive loop: 6: for subiteration=1 to 2 θ + 1 do

Figure 4

 4 Figure4shows a DAG for an iteration with a maximum temporal level θ = 2, three temporal levels of cells, 4 subiterations 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.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.

Fig. 5 :

 5 Fig. 5: Mesh for the take-off blast wave computation.

Fig. 6 :

 6 Fig. 6: Take-off blast wave computation.

Fig. 7 :

 7 Fig. 7: Impact of Priority Strategy with varying number of CEs. First, let us consider the computation with 128 CEs to highlight what happens when we use or not the priority strategy. Figure 8 is a Gantt Diagram where each horizontal bar represents a worker. Idle time is always in red, the overhead is in yellow.The first step of the computation is in green (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 priority they

  Fig. 8: 128 CE computation without or with priorities.

Fig. 9 :

 9 Fig. 9: Evolution of ready tasks (Solver step, 128 CEs).

Fig. 10 :

 10 Fig. 10: Time performance with different configurations with parallel workers.

(a) 8 Fig. 11 :

 811 Fig. 11: Time spent in the solver with different configurations (X-axis gives the number of nodes).

TABLE I :

 I Cell distribution and associated computation.

TABLE II :

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

	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

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


Thomas Simon

ON THE HAUSDORFF DIMENSION OF REGULAR POINTS OF INVISCID BURGERS EQUATION WITH STABLE INITIAL DATA

Keywords: 2000 Mathematics Subject Classification. 35Q53, 35R60, 60F99, 60G52, 60J30 Burgers equation, Hausdorff dimension, integrated stable process, lower tail probabilities, shock structure

. Along the way, we contradict a recent conjecture of Z. Shi about the lower tails of integrated stable processes.

Introduction and statement of the results

Since the seminal paper [START_REF] Ya | Statistics of shocks in solutions of inviscid Burgers equation[END_REF], statistical properties of the Burgers equation where ψ 0 is the initial potential given by u 0 = -∂ x ψ 0 . The inviscid Burgers equation is a simplified version of (1.1) where the viscosity parameter ν = 0, and its so-called Hopf-Cole solution is obtained from (1.2) in letting ν → 0. By Laplace approximation, it takes a particularly nice form:

(1.3) ψ(t, x) = sup a∈R ψ 0 (a) - (x -a) 2 2t ,
which is well-defined provided the initial potential satisfies ψ 0 (a) = o(a 2 ) when |a| → +∞.

We refer e.g. to the monograph [START_REF] Woyczynski | Burgers-KPZ turbulence[END_REF] for the above facts, and much more, concerning Burgers equation.

In this paper, we are interested in the inviscid Burgers equation whose initial data is a two-sided Lévy α-stable process. More precisely, we suppose that the initial condition X is defined followingly:

(1.4)

X x = Z x if x ≥ 0 -Z ′ -x if
x ≤ 0 where Z = {Z x , x ≥ 0} is a Lévy α-stable process and Z ′ an independent copy of Z. Referring to Chapter VIII in [START_REF] Bertoin | Lévy Processes[END_REF] for more details, let us recall that Z is a real process with stationary and independent increments, which is (1/α)-self-similar:

(1.5) {Z kx , x ≥ 0} d = {k 1/α X x , x ≥ 0}
for all k > 0. This property forces the stability index α to be in (0, 2], and Z is Brownian motion (up to a scaling parameter) in the case α = 2. Its Lévy-Khintchine exponent Ψ(λ) = log E[e iλZ 1 ] is given by

Ψ(λ) = κ|λ| α (1 -iβsgn(λ) tan(πα/2)), λ ∈ R,
where κ > 0 is the scaling parameter and β ∈ [-1, 1] is the skewness parameter. Without loss of generality, in the following we will suppose κ ≡ 1. The positivity parameter ρ = P [Z 1 > 0] takes its values in [1-1/α, 1/α] when α > 1 and in [0, 1] when α < 1. When α > 1, the value ρ = 1 -1/α corresponds to the spectrally positive situation (β = 1 and Z has no negative jumps) and the value ρ = 1/α to the spectrally negative situation (β = -1 and Z has no positive jumps). When α < 1 and ρ = 0 (resp. α < 1 and ρ = 1), Z is the negative of a subordinator and has no positive jumps (resp. a subordinator and has no negative jumps). When α = 1, then ρ ∈ (0, 1) and Z has jumps in both directions.

The law of the iterated logarithm for Z at infinity -see Theorem VIII.5 in [START_REF] Bertoin | Lévy Processes[END_REF] -entails

lim sup x→+∞ x -κ Z x = 0 or + ∞ according as κ > 1/α or κ ≤ 1/α.
Hence, since our initial potential is given up to some meaningless additive constant by

ψ 0 (x) = x 0 X t dt, x ∈ R,
the growth condition ψ 0 (a) = o(a 2 ) at infinity assigns the restriction

α ∈ (1, 2]
on the stability parameter, which will be supposed henceforth unless explicitly stated.

Derivating (1.3) with respect to x yields readily the following formula for the Hopf-Cole solution of (1.1):

u(t, x) =
xa(t, x) t where a(t, x) is the largest point attaining the maximum in (1.3), in other words:

a(t, x) = max s ∈ R, C ′ t (s) ≤ xt -1
where C ′ t stands for the right-derivative of C t , which is the convex hull of the function

x → x 0 (X u + ut -1 )du.
The so-called Lagrangian regular points of (1.1) are the points where the above function coincides with its convex hull. Notice that this time-dependent set L t can be described in terms of the function a(t, x):

L t = {a(x, t), x ∈ R and a(x-, t) = a(x, t)} .
From the physical point of view, L t is the set of particles which have not been shocked up to time t by the turbulence governed by Equation (1.1) -see [START_REF] Woyczynski | Burgers-KPZ turbulence[END_REF] -and for this reason there has been some interest over the years in describing the structure of the sets L t , especially in studying their fractal properties. In [START_REF] Janicki | Hausdorff dimension of regular points in stochastic Burgers flows with Lévy α-stable initial data[END_REF], the authors raised the following Conjecture (Janicki and Woyczynski). For every t > 0, one has

Dim H L t = 1/α a.s.
It is easy to see that the Hausdorff dimension of L t does not, indeed, depend on t in this model. Namely, from the self-similarity of Z one can show -see [START_REF] Janicki | Hausdorff dimension of regular points in stochastic Burgers flows with Lévy α-stable initial data[END_REF] p. 285 -that u(t, x)

d = t 1/(α-1) u(1, xt α/(α-1) ),
which entails that a.s. Dim H L t = Dim H L 1 for every t > 0. In the following, we will be therefore interested in the set L 1 only, which we denote by L for the sake of simplicity. Notice that the above conjecture had been previously solved (without complete proof) by Sinai [START_REF] Ya | Statistics of shocks in solutions of inviscid Burgers equation[END_REF] in the Brownian case α = 2, and that Bertoin [START_REF] Bertoin | The inviscid Burgers equation with Brownian initial velocity[END_REF] showed then rigourously the result in the general case α ∈ (1, 2] with no positive jumps, as a consequence of the remarkable fact that the process x → a(x, 1) is a subordinator close to the(1/α)-stable onesee Theorem 2 in [START_REF] Bertoin | The inviscid Burgers equation with Brownian initial velocity[END_REF]. Nevertheless, this latter property is no more true in the non spectrally negative framework -see the conclusion of [START_REF] Bertoin | The inviscid Burgers equation with Brownian initial velocity[END_REF], and the structural study of (1.1) seems to require an entirely different methodology when there are positive jumps. To this end, let us also mention an attempt made in [START_REF] Carraro | Equation de Burgers avec conditions initiales à accroissements indépendants et homogènes[END_REF] with the concept of statistical solution (which is different from the Hopf-Cole solution).

In this paper we show that Janicki and Woyczynski's conjecture is false in general when there are positive jumps:

Theorem A. Suppose that Z has positive jumps. Then there exists α 0 > 1 such that

Dim H L < 1/α a.s.
for every α ∈ (1, α 0 ).

Our main argument comes from a recent paper by Molchan and Kholkhov [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF], who were interested in the sets L t of Lagrangian regular points of (1.1) when the initial data is a two-sided fractional Brownian motion W . In [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF], they proved that an a.s. upper bound on Dim H L t -which is also independent of t > 0 by the self-similarity of W -follows from a lower bound on the exponent κ appearing in the estimate

P t 0 Ŵs ds < ε + t 2 , ∀ t ∈ [-1, 1] ≤ ε κ , ε → 0, where Ŵ = -W d = W
is the dual process of W . We remark that this argument transfers to the Lévy stable case without much difficulty, and is actually even simpler because of the independence and stationarity of the increments of X. This will be done in Section 3. Before this, we will have to prove a crucial estimate on the first-passage time of the integrated stable process. Let us fix first some more notation, introducing

A t = t 0 Z s ds, t ≥ 0,
the integral of a Lévy α-stable process Z -with the same notations as above, but this time for every α ∈ (0, 2] -and T = inf{t > 0, A t = 1} the first-passage time of A across 1.

Theorem B. Suppose that Z has negative jumps and is not the negative of a subordinator. With the above notations, for every α 0 > 0 there exists a constant κ > 0 depending only on ρ such that

lim inf t→∞ t κ P[T > t] = 0 for every α ∈ [α 0 , 2].
This result, which will be proved in Section 2, is interesting in its own right because it contradicts another conjecture whose solution had been announced (with a hidden error) by the author in [START_REF] Simon | The lower tail problem for the area of a symmetric stable process[END_REF], and whose statement was the following: Conjecture (Shi). Suppose that Z is symmetric and that α > 1. Then with the above notations,

(1.6) P[T > t] = t -(α-1)/2α+o(1) , t → ∞.
Indeed, since in Theorem B the exponent κ is positive independent of α, we see that the above conjecture is contradicted when α is close enough to 1, as in Theorem A. Notice that in a recent paper [START_REF] Simon | The lower tail problem for homogeneous functionals of stable processes with no negative jumps[END_REF] the author proved (1.6) in the case where Z has no negative jumps and α > 1. The idea -which had been originally given by Z. Shi -consisted in time-changing the process A through τ the inverse local time of Z and considering the fluctuations of the Lévy stable process A τ . We have been thinking for a long time that this method would be also successful when Z has negative jumps, because A τ is a Lévy symmetric (α-1)/(α+1)-stable process whatever the value of ρ should be -see Lemma 1 in [START_REF] Simon | The lower tail problem for homogeneous functionals of stable processes with no negative jumps[END_REF], which allows in particular to obtain a general upper bound

P[T > t] ≤ K t -(α-1)/2α , t → ∞ for every value of ρ ∈ [1 -1/α, 1/α]
and some finite constant K -see Theorem A in [START_REF] Simon | The lower tail problem for homogeneous functionals of stable processes with no negative jumps[END_REF].

It now appears that this subordination method is too crude when there are negative jumps, and yields only an upper bound which is not optimal, at least when α is close to 1. In this paper, we will obtain a uniform upper bound by discretization through exponential timechange combined with FKG-type inequalities, all of which we learned from Caravenna and Deuschel in the genesis of their paper [START_REF] Caravenna | Pinning and wetting transition for (1+1)-dimensional fields with Laplacian interaction[END_REF]. Let us stress that these arguments are also by far non optimal. Nevertheless they are quite robust and, since they involve eventually only fixed upper tails of Z and A which can be bounded independently of α, this method makes it at least possible to contradict both Janicki-Woyczynski and Shi's conjectures when α is close to 1. Actually, we believe that these conjectures are false for all values of α except when there are no positive jumps (for the first) or no negative jumps (for the second), and in the fourth and last section of the paper we will state two other conjectures for the values of Dim H L and the critical exponent in (1.6), in a general non completely asymmetric framework.

Proof of Theorem B

We begin with the case α > 1, and we will actually obtain a slightly stronger result which is crucial for Theorem A. Fix ρ = P[Z 1 > 0] > 0 once and for all, and set γ = (α -1)/α > 0. We will show that there exists κ > 0 independent of α such that

(2.1) lim inf t→+∞ t κ P A s < 1 + s + t -γ s 2 , ∀ s ≤ t = 0,
which readily entails Theorem B by comparison. We first define an exponential subsequence of times, in considering the events

A n = A 2 m < 1 + 2 m + 4 m-nγ , m = 0 . . . 2n , n ≥ 0.
Clearly, it is enough to prove that there exists κ > 0 independent of α such that P[A n ] ≤ e -κn for n sufficiently large. We will obtain slightly more, in showing that (2.2)

P[B n ] ≤ e -κn
for some κ > 0 independent of α and n large enough, where

B n = A 2 m < 1 + 2 m + 4 m-nγ , m = 0 . . . n , n ≥ 0.
To do so, consider the events

C k = {Z 2 k > -2 k/α , A 2 k > -2 k(1+1/α)
} for every k ≥ 0, and the family of random times defined recursively by

σ 0 = 0 and σ n = inf{k > σ n-1 / C k occurs}, n ≥ 1.
If {F t , t ≥ 0} stands for the completed σ-field generated by {Z s , s ≤ t}, then 2 σn is a F t -stopping time for every n ≥ 0. Denoting henceforth by [x] the integer part of any real number x, we now state a crucial lemma which is mainly borrowed from [START_REF] Caravenna | Pinning and wetting transition for (1+1)-dimensional fields with Laplacian interaction[END_REF]:

Lemma. For every α 0 > 0, there exists δ, c > 0 independent of α ∈ [α 0 , 2] such that P σ [δn] < n | B n ≥ c
for all n sufficiently large.

Taking this lemma for granted and setting K = c -1 < +∞ we see that

P [B n ] ≤ KP σ [δn] < n, B n ≤ KP A 2 m < 1 + 2 m + 4 m-nγ ∀ m = σ 1 + 1, . . . , σ [δn] + 1, σ [δn] < n .
Introducing the notations

G [δn] = F σ [δn] and D [δn] = A 2 m < 1 + 2 m + 4 m-nγ ∀ m = σ 1 + 1, . . . , σ [δn] + 1
for every n ≥ 0, the strong Markov property at time 2 σ [δn] yields

P D [δn] , σ [δn] < n ≤ P 1 D [δn]-1 ∩{σ [δn] <n} P A 2 σ [δn] +1 < 1 + 2 σ [δn] +1 + 4 σ [δn] +1-nγ | G [δn] ≤ P 1 D [δn]-1 ∩{σ [δn] <n} P A 2 σ [δn] +1 < 1 + 2 σ [δn] +1 + 2 2+2σ [δn] /α | G [δn]
where in the second line we used the event {σ [δn] < n}. On the other hand, again by the strong Markov property, conditionally on G [δn] we can decompose

A 2 σ [δn] +1 = A 2 σ [δn] + 2 σ [δn] Z 2 σ [δn] + Â2 σ [δn]
where  is a copy of A independent of G [δn] . By the definition of σ [δn] , this yields

P[A 2 σ [δn] +1 < 1 + 2 σ [δn] +1 + 2 2+2σ [δn] /α | G [δn] ] ≤ P[ Â2 σ [δn] < 2 2+σ [δn] (1+1/α) + 2 2+2σ [δn] /α | G [δn] ] ≤ P[ Â2 σ [δn] < 2 3+σ [δn] (1+1/α) | G [δn] ] ≤ P[A 1 < 8],
where we used the fact that α > 1 in the second line, and the (1 + 1/α)-self-similarity of  independent of G [δn] in the third. But from Proposition 3.4.1 in [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF], we know that A 1 is a real α-stable random variable whose Lévy-Khintchine exponent Φ(λ) =log E[e iλA 1 ] is given by

Φ(λ) = (α + 1) -1 |λ| α (1 -iβsgn(λ) tan(πα/2)), λ ∈ R,
where β is the skewness parameter of Z 1 . In particular, its positivity parameter is P[A 1 > 0] = ρ > 0 and its scaling parameter belongs to [1/3, 1/2) independent of α. This clearly entails -see Section 1.6 in [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF] -that there exists κ < 1 depending only on ρ such that

P[A 1 ≥ 8] ≥ (1 -κ) for every α > 1.
Let me also stress that this argument breaks down when there are no negative jumps, since then P[A 1 ≥ 8] → 0 when α → 1. We finally obtain

P D [δn] , σ [δn] < n ≤ κP D [δn]-1 , σ [δn] < n ≤ κP D [δn]-1 , σ [δn]-1 < n -1 ≤ κ [δn]
by an induction argument, and since above the event {σ [δn] < n} is only used to obtain the upper bound 4 σ [δn] +1-nγ ≤ 2 2+2σ [δn] /α . Now because δ, κ are independent of α, this yields (2.2) as desired, and completes the proof of Theorem B when α > 1. The case α ≤ 1 can be handled exactly in the same way in working on the events

{A 2 m < 1 + 2 m , m = 0 . . . n} , n ≥ 0,
and we leave the details to the reader. However, for the sake of completeness and since our arguments are partly different from [START_REF] Caravenna | Pinning and wetting transition for (1+1)-dimensional fields with Laplacian interaction[END_REF], we will give the Proof of the lemma. Set P n [.] = P[ . |B n ] for concision. By the definition of σ [δn] , we have for every δ > 0

P n [σ [δn] < n] = 1 -P n n k=1 1 C c k > δ ′ n ≥ 1 - 1 δ ′ n n k=1 P[C c k ] = 1 δ ′ n n k=1 P n [C k ] - δ δ ′
where we set δ ′ = 1δ for conciseness. Hence, it suffices to prove the existence of c > 0 independent of α and k such that

P n [C k ] ≥ c
for every k ∈ [1, n] and every n sufficiently large. To do so, we will first consider the events

D k,n = Z t < 2 -(n-k)/α , ∀ t ∈ [0, 2 n+k ] ⊆ B n
and prove the intuitively obvious inequalities

(2.3) P[C k | B n ] ≥ P[C k | D k,n ]
for every n ≥ 1 and every k ∈ [1, n], with the help of a discretization of Z and a FKG argument. Fixing k and n, set

Z N t = i<tN Li+1 N -L i N and A N t = 1 N i<tN j≤i Lj+1 N -L j N for every t ≥ 0, N ≥ 1.
P[C N k | B N n ] ≥ P[C N k | D N k,n ]
for every N ≥ 1. The key-point is that the probabilities of the events B N n , C N k , D N k,n and their intersections depend only on the joint law P N of the increments {(Z i N -Z i-1 N ), i = 1, . . . , N2 n+k }, which are stationary and independent. The measure P N has the density

(2.5) f (x) = N 2 n+k i=1 g(x i )
with respect to the Lebesgue measure on R N 2 n+k , where g is the density of the variable Z 1 N , and the density of the conditional measure P N n [.] = P N [ . | BN n ] with respect to the Lebesgue measure on R N 2 n+k is hence given by

h(x) = 1 BN n f (x) P[B N n ] where BN
n is the image of B N n under P N . Besides, it is obvious from the definition of B N n that the function 1 BN n (x) is decreasing, in the sense that if x i ≥ y i for all i = 1 . . . N2 n+k , then 1 BN n (x) ≤ 1 BN n (y). From (2.5) and this monotonicity property, we deduce that h satisfies Holley's criterion:

h(x)h(y) ≥ h(x ∨ y)h(x ∧ y), x, y ∈ R N 2 n+k
with the notations (x∨y) i = x i ∨y i and (x∧y) i = x i ∧y i for all i = 1 . . . N2 n+k . By Corollary 12 in [START_REF] Holley | Remarks on the FKG inequalities[END_REF], this entails that the measure P N n satisfies FKG inequality in the sense that n . This finally entails

P N n [C ∩ D] ≥ P N n [C]P N n [D]
P N n [ ĈN k ∩ DN n ] ≤ P N n [ ĈN k ]P N n [ DN n ]
and, after some elementary transformations using the inclusion D N k,n ⊆ B N n , the desired inequality (2.4) for every k = 1, . . . , n and n, N ≥ 1. Hence, from (2.3), we now need to show that there exists c > 0 independent of α such that

P[C k | D k,n ] ≥ c
foe every k ∈ [1, n] and n sufficiently large. A scaling argument yields first

P[C k | D k,n ] = P[Z 1 < 1, A 1 < 1 | {Z t > -2 -n , ∀t ∈ [0, 2 n ]}]
where P stands for the law of the dual process Ẑ = -Z. By Chaumont's results -see Remark 1 and Theorem 6 in [START_REF] Chaumont | Conditionings and path decomposition for Lévy processes[END_REF], the conditional law on the right-hand side converges to P↑ which is the law of Ẑ conditioned to stay positive. Hence, for n sufficiently large, one has

P[C k | D k,n ] ≥ 1 2 P↑ [Z 1 < 1, A 1 < 1] ≥ 1 2 P↑ [S 1 < 1]
with the notation S 1 = sup{Z t , t ≤ 1}. It is now intuitively obvious that the right-hand side can be bounded from below by a positive constant on every closed interval [α 0 , 2] ⊂ (0, 2].

Let us give a rigorous argument. Using Theorem 1 in [START_REF] Chaumont | Excursion normalisée, méandre et pont pour les processus de Lévy stables[END_REF], one gets

P↑ [S 1 < 1] ≥ c Ê(me) [Z 2 1 1 {S 1 <1}
] for some constant c depending only on ρ, where P(me) denotes the law of the meander associated to Z. The pathwise representation of the meander given in [START_REF] Bertoin | Lévy Processes[END_REF], Proposition VIII. [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF] yields

Ê(me) [Z 2 1 1 {S 1 <1} ] ≥ c ′ Ê[(Z 1 -I 1 ) 2 1 {S 1 -I 1 <1}
] for some constant c ′ depending only on ρ, with the notation I 1 = inf{Z t , t ≤ 1}. Finally, the Lévy-Khintchine formula and Theorem 1 in [START_REF] Kadankova | On the joint distribution of the supremum, infimum, and the value of a semicontinuous process with independent increments[END_REF] entail that the probability on the righthand side is a continuous function of the parameter α. Since it is obviously positive for every α ∈ (0, 2], it is bounded from below by a positive constant on every closed interval [α 0 , 2] ⊂ (0, 2], which completes the proof.

Remark. The proof of the above lemma would be slightly simpler if we could prove that the measures P n themselves satisfy FKG inequality, instead of considering their discretizations P N n . From Example 2.3.6 and Theorem 4.6.1 in [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF], and by right-continuity of the sample paths of Z, we know that the unconditioned measure P satisfies FKG, in the sense for every time-horizont T > 0 and every bounded measurable increasing functionals F, G :

D T → R + , (2.6) E[F (Z t , t ≤ T )G(Z t , t ≤ T )] ≥ E[F (Z t , t ≤ T )]E[G(Z t , t ≤ T )]
(here we set D T for the Skorokhod space of càdlàg functions from [0, T ] to R, and we say that a functional F :

D T → R + is increasing if for every x, y ∈ D T , x t ≥ y t ∀ t ≤ T =⇒ F (x t , t ≤ T ) ≥ F (y t , t ≤ T ).) Since B n is a monotonous event, our desired conditioned version of (2.6) (2.7) E n [F (Z t , t ≤ T )G(Z t , t ≤ T )] ≥ E n [F (Z t , t ≤ T )]E n [G(Z t , t ≤ T )]
would be fulfilled if P satisfied the strong FKG inequality. However, there exist some pathmeasures which are FKG but not strong FKG -we learned this from J.-D. Deuschel, and we could not prove (2.7) directly, even in the Brownian case.

Proof of Theorem A

Recall that we are interested in the Hausdorff dimension of the random set

L = {a(x), x ∈ R and a(x-) = a(x)}
where a(x) := max{s ≥ 0, C ′ (s) ≤ x} and C ′ is the right-derivative of C, the convex hull of the function

x → x 0 (X u + u) du, x ∈ R.
Recall also that X is a two-sided α-stable Lévy process (α > 1) with positive jumps as defined in (1.4), and fix its positivity parameter ρ < 1/α once and for all. Notice that by definition X does not jump negatively at time l whenever l ∈ L: X l ≥ X l- a.s. On the other hand, it is well possible that l ∈ L is a "conical point" in the sense that X l > X l-. However, if we define

L = a(x), x ∈ R, a(x-) = a(x) and X a(x) = X a(x)-,
we see from the fact that the set of points of discontinuity of X is a.s. countable that

Dim H L = Dim H L a.s.
The key-point -which was first noticed by Sinai [START_REF] Ya | Statistics of shocks in solutions of inviscid Burgers equation[END_REF] in the Brownian case -is that a.s.

L ⊆ L := a ∈ R / x 0 (X u + u) du ≥ a 0 (X u + u) du + (x -a)(X a + a), ∀ x ∈ R ,
so that we only need to get an upper bound on Dim H L. To do so, we will use (2.1) together with the same arguments as Molchan and Kholkhov [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF]. First, a slight modification of Lemma 1 in [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF] shows that Theorem A will be proved as soon as there exists a constant κ > 0 independent of α and a subsequence δ n → 0 such that

P L ∩ (x -δ n , x + δ n ) = ∅ ≤ δ κ n , n → ∞
uniformly in x ∈ R. Indeed, reasoning exactly as In Lemma 1 in [START_REF] Molchan | Small values of the maximum for the integral of fractional Brownian motion[END_REF] entails then

Dim H L ≤ 1 -κ a.s.
and completes the proof of the theorem with α 0 = 1/(1κ). Second, we remark that by linearity of the integral and by independence and stationarity of the increments of the process u → X u + u, the random sets

L ∩ (x -δ n , x + δ n ), x ∈ R
have all the same law. Hence, we need to show that there exists a constant κ > 0 independent of α such that (3.1) lim inf δ→0 δ -κ P Lδ = ∅ = 0, where for every δ > 0 we wrote

Lδ = |a| < δ / x 0 (X u + u) du ≥ a 0 (X u + u) du + (x -a)(X a + a), ∀ x ∈ R .
For every δ > 0, set F δ for the completed filtration generated by {X u , |u| ≤ δ}, and consider the F δ -measurable random variables

M δ = sup |x|≤δ |X x + x| and N δ = sup |x|≤δ x 0 (X u + u) du .
If X = -X and Ẑ = -Z denote the dual processes of X and Z respectively, we see that a.s.

Lδ = ∅ ⊆ x 0 ( Xu -u) du ≤ N δ + (δ + |x|)M δ ∀ x ∈ R ⊆ A δ ∩ A ′ δ
with the notations

A δ = x δ (( Ẑu -Ẑδ ) -(u -δ)) du ≤ 2N δ + 2(δ + x)M δ ∀ x ≥ δ and A ′ δ = -x -δ ((-Ẑ′ u + Ẑ′ δ ) -(u -δ)) du ≤ 2N δ + 2(δ + |x|)M δ ∀ x ≥ δ ,
where Ẑ′ is an independent copy of Ẑ. By the Markov property, the events A δ and A ′ δ are independent and identically distributed conditionally on F δ , so that (3.2)

P Lδ = ∅ = P P Lδ = ∅ | F δ ≤ E [P [A δ ∩ A ′ δ | F δ ]] = E P [A δ | F δ ] 2 .
Besides, again by the Markov property, we can write

P [A δ | F δ ] = P t 0 L u du ≤ 2n + 2tm + t 2 /2, ∀ t ≥ 0 n=N δ ,m=M δ
where L is a copy of Ẑ independent of F δ . Now since α > 1 and by the scaling property of X, it is immediate to see that a.s. N δ ≤ δ 1+1/α (1 + N) and M δ ≤ δ 1/α (1 + M) as soon as δ < 1, where M, N are F δ -measurable and such that

M d = sup |x|≤1 |X x | and N d = sup |x|≤1 x 0 X u du .
Setting ε = (2δ) (α+1)/2α and R = max{(1 + N), (1 + M) α+1 } for conciseness, we can rewrite

P [A δ | F δ ] = P t 0 L u du ≤ ε 2 r + (ε 2 r) 1/(α+1) t + t 2 /2, ∀ t ∈ [0, 1] r=R •
Returning to (3.2), we obtain

P Lδ = ∅ ≤ P R ≥ ε -1 + P t 0 L u du ≤ ε + ε 1/(α+1) t + t 2 /2, ∀ t ∈ [0, 1] 2 ≤ (cδ) 1/4 + P t 0 L u du ≤ ε + ε 1/(α+1) t + t 2 , ∀ t ∈ [0, 1] 2
for some c > 0, where in the third line we used well-known estimates on the upper tails of supremum of stable processes -see e.g. Theorem 10.5.1. in [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF]. A scaling argument yields

P t 0 L u du ≤ ε + ε 1/(α+1) t + t 2 , ∀ t ∈ [0, 1] = P t 0 L u du ≤ 1 + t + n -γ t 2 , ∀ t ∈ [0, n]
where n = ε -α/(α+1) = 1/ √ 2δ. Finally, since L d = -Z has negative jumps, we see from (2.1) that there exists κ > 0 depending only on ρ such that lim inf δ→0 δ -κ P Lδ = ∅ = 0, which is (3.1) and completes the proof.

Two conjectures

Let us start by a classical result of Bingham concerning the asymptotics of the ruin probabilities related to the stable process Z: if we set S = inf{t > 0, Z t > 1}, then there exists a constant c ∈ (0, ∞) such that (4.1)

P[S > t] ∼ ct -ρ , t → +∞
as soon as |Z| is not a subordinator -see Proposition VIII. 2 in [START_REF] Bertoin | Lévy Processes[END_REF]. Notice in passing that the constant c equals α p Z 1 (0) when Z has no positive jumps -this is a consequence of Skorokhod's formula written e.g. in [START_REF] Bingham | Fluctuation theory in continuous time[END_REF] p. 749, and that in the other cases it can be given (non explicitly) in terms of the excursion measure associated to the reflected process -see Lemma 1 in [START_REF] Chaumont | Excursion normalisée, méandre et pont pour les processus de Lévy stables[END_REF]. Considering now T = inf{t > 0, A t = 1} the first-passage time of the integral of Z across 1, from (4.1) it is tantalizing to state the Conjecture C. Suppose that |Z| is not a subordinator. Then

P[T > t] = t -ρ/2+o(1) , t → +∞.
The Brownian case α = 2 and ρ = 1/2 had been obtained in [START_REF] Goldman | On the first passage of the integrated Wiener process[END_REF] after expanding a closed formula of McKean concerning the distribution of T , which yields actually a more accurate estimate like (4.1) with an explicit constant for c -see Proposition 2 therein. The case with no negative jumps α > 1 and ρ = (α -1)/α was proved recently in [START_REF] Simon | The lower tail problem for homogeneous functionals of stable processes with no negative jumps[END_REF], with a good control on o(1) allowing to show that E[T ρ/2 ] = +∞. The above conjecture is also motivated by the aforementioned fact that Z and A have the same positivity parameter ρ, a quantity which should typically play a rôle in the distribution of T . On the other hand, the process A is smoother than Z, so that the upper tails of the distribution of T should be heavier than those of S. We propose the value ρ/2 for the critical exponent, since it is in accordance with the spectrally positive case. Conjecture C will be the matter of further research [START_REF] Devulder | The lower tail problem for integrated stable processes[END_REF].

Suppose now α > 1 and consider the drifted stable process Z c t = Z t + ct, t ≥ 0, for some c = 0. From Lemma VI.21 in [START_REF] Bertoin | Lévy Processes[END_REF] and explicit estimates on the renewal function of Z c , one can show that

P[Z c t < ε, ∀ t ∈ [0, 1]] ≍ ε ρα , ε → 0
for every c ∈ R. Notice also that the latter estimate is false when α = 1 by Bingham's result, since Z c is then a strictly stable process whose positivity parameter depends on c. We think that the estimate is also untrue when α < 1, but we got stuck in proving this. Setting A c

for the integral of Z c , we believe from the above fact that when α > 1, (4.2) P[A c t < ε, ∀ t ∈ [0, 1]] = ε ρα/(α+1)+o (1) , ε → 0 for every c ∈ R. By self-similarity, notice that this estimate is the same as Conjecture C when c = 0. When c < 0, it is particularly relevant for the inviscid Burgers equation whose initial data is the two-sided dual process Ẑ. Indeed, reasoning as in the proof of Theorem A, one can show that the upper bound in (4.2) entails Dim H L ≤ ρ a.s. where ρ = 1ρ is the positivity parameter of Ẑ. Considering now (1.1) where the initial data is the two-sided dual process Z, from Conjecture C, the above considerations and Bertoin's result [START_REF] Bertoin | The inviscid Burgers equation with Brownian initial velocity[END_REF], one is tempted to state the Conjecture D. With the above notations, for every α > 1 and every ρ ∈ [1 -1/α, 1/α] one has Dim H L = ρ a.s.

From the present paper, we are convinced that optimal lower tail estimates for the integral of Ẑ should provide the key-argument to obtain the upper bound in Conjecture D. The lower bound seems more delicate because of the positive jumps of Ẑ which prevent from using Handa's criterion [START_REF] Handa | A remark on shocks in inviscid Burgers' turbulence[END_REF]. On the other hand, it is immediate to see from its definition that the set L t contains the points of global increase of the drifted process Z 1/t for every t > 0, and one may wonder if the methods developped by Marsalle [START_REF] Marsalle | Hausdorff measures and capacities for increase times of stable processes[END_REF] to determine the Hausdorff dimension of the points of local increase of the non-drifted process Z, could not be useful.

(1. 1 ) 2 R exp 1 2ν ψ 0

 120 ∂ t u + u∂ x u = ν∂ xx u, ν > 0 with initial condition u 0 (x) := u(0, x) = X x where {X x , x ∈ R} is a given random process, have given rise to intensive research. Eventhough (1.1) is a much simplified version of the Navier-Stokes equation, it is still relevant in physics as a model equation for e.g. shock waves in hydrodynamics, in order to describe turbulence phenomena not covered by the classical Korteweg-de Vries equation. From the mathematical point of view, a nice feature of (1.1) is the possibility to solve it explicitly through the change of variable u = -∂ x ψ, which is known as the Hopf-Cole substitution: one has (1.2) ψ(t, x) = 2ν log (4πνt) -1/(a) -(xa) 2 2t da

whenever 1 C

 1 and 1 D are increasing functions. On the other hand, the function 1 ĈN k (x) is increasing and the function 1 DN n (x) is decreasing, with the same notation as above for ĈN k and DN

  n+k ] for the Skorokhod topology when N → ∞. Hence setting B N n , C N k , D N k,n for the events B n , C k , D k,n with (Z, A) replaced by (Z N , A N ), by weak convergence and right continuity it suffices to show

	(2.4)

It is well-known -see e.g. Chapter 3 in [3] -that the bivariate process {(Z N t , A N t ), t ∈ [0, 2 n+k ]} converges in law towards (Z t , A t ), t ∈ [0, 2

Acknowledgement. I am very grateful to F. Caravenna and J.-D. Deuschel for having let me experience, live in Berlin, the ongoing of their paper [5].


Laura Pometcu
email: laura.pometcu@univ-rennes1.fr.

Ala Sharaiha

Ratiba Benzerga

Razvan D Tamas

Philippe Pouliguen

) A Sharaiha

Method for material characterization in a non-anechoic environment

I. INTRODUCTION

Numerous methods of obtaining the dielectric properties of materials at microwave frequencies have been proposed in the past. The standard measurement methods are the co-axial probe, transmission line, cavity, free space methods, and so on [START_REF] Chen | Microwave Electronics: Measurement and Materials Characterization[END_REF][START_REF] Fiorillo | [END_REF][3][START_REF]Basics of measuring the dielectric properties of materials[END_REF][START_REF] Bussey | Proceeding of the IEEE[END_REF][START_REF] Fenner | [END_REF][7] . Nonetheless, these methods have their advantages and disadvantages depending on the application. Waveguide methods are popular, where the sample is precisely machined to fit inside the waveguide. Both rectangular and coaxial waveguides are used, rectangular samples are easier to produce than coaxial ones, however they can only be used over a limited frequency range. A coaxial waveguide allows for extremely wide-band measurements but the sample preparation is difficult. Free-space techniques [8][START_REF]Solutions for Measuring Permittivity and Permeability with LCR Meters and 14 Impedance Analyzers[END_REF][START_REF] Schultz | [END_REF][START_REF] Zivkovic | [END_REF] circumvent the problem of the sample fit precision and are nondestructive to the sample. A major disadvantage of the free-space method is that it is not done in an enclosed space, and therefore a greater sample surface needs to be fabricated to avoid diffraction effects around the sample edges. Some reductions of the required sample size for free space microwave measurements have been demonstrated by incorporating spot-focusing horn antennas, 8,12 but this technique is limited in bandwidth due to the focusing nature of the lenses.

As opposed to the classical techniques that use an anechoic chamber, measuring the reflection coefficient for a material sample in a non-anechoic environment should take into consideration a multipath transmission scenario (including reflection scattering and diffraction on the walls and on other objects encountered) 13 .

A post processing after the measurement is necessary in order to extract accurately the reflection coefficient of the material sample under test. This is called a "de-embedding" technique [START_REF] Wee | International Microwave and Optoelectronics Conference[END_REF] . In literature, most of the characterization methods require a sample of a size much greater than the antenna aperture to be used, as most of the wave front surface should be incident on the material sample 8,[START_REF] Ghodgaonkar | Instrumentation and Measurement[END_REF] .

For characterizing a material sample with small dimensions and at low frequencies (where the size of the antennas is very important) the classical "de-embedding" method cannot be applied 8,[START_REF] Hanson | Instrumentation and Measurement[END_REF] .

In this paper we propose a method for measuring the reflection coefficient of materials that does not require large material samples and is done in a non-anechoic environment. It is based on combining the "de-embedding" and the average value techniques for several sample distances, resulting in the elimination of unwanted multiple diffractions and reflections for the environment. The material's properties can then be deduced from the reflection coefficient. This coefficient is extracted from the transmission between two antennas through a bistatic reflection method 17 .

II. PROPOSED APPROACH AND MEASUREMENT SETUP

A schematic diagram of the bistatic microwave measurement system is given in Fig. 1. The transmitting and receiving antennas are focusing broadband horn antennas QWH-SL-2-18-N-HG-R 18 .

They are placed side by side and oriented towards a Material Under Test (MUT). A vector network analyzer is used for this type of measurements. The size of the MUT depends on the desired frequency range. At low frequencies we would need large samples and smaller distances between the antennas and the samples in order to eliminate the unwanted reflections. The calibration has been done only for the coaxial cables using an OSLT (open, short circuit, load and through line) calibration kit. Since the antennas are close one to the other the angle between the material sample and the two antennas Ф is close to 0 0 . Thus, we can do the approximation that the transmission coefficient S 21 measured with the VNA corresponds to the reflection coefficient of the MUT but at the reference point X. We need the reflection coefficient at the surface of the material, in Fig. 1 referred to as Γ. The parameter S 21 can be obtained by calculating the ratio between the reflected electric field over the incident electric field on the MUT of a planar wave at a distance R from the MUT in far field [START_REF] Balanis | Antenna Theory Analysis and Design Third Edition[END_REF][START_REF] Balanis | Advanced Engineering Electromagnetics[END_REF] , such as: (1) where A is the mutual coupling between the two antennas and the medium itself, B the amplitude of the incident wave [START_REF] Xu | Silicon Monolithic Integrated Circuits in RF Systems[END_REF] , k 0 the wave number, R the distance between the antennas and the sample and Γ MUT is the reflection coefficient of the Material. The interest is focused on the reflection coefficient of the material, Γ MUT that will be extracted using equation (4).

A measurement with a metallic plate in front of the antennas that replaces the material sample and another measurement with the antennas alone (without the material sample or the metallic plate) are performed. These two other measurements are necessary in order to establish a reference. For these two cases, the transmission coefficients are:

(2) and ,

(3) respectively.

From equations ( 1) -( 3) we can extract the exact reflection coefficient of the MUT using the following formula: (4) This formula can be considered as the "de-embedding" process.

The characterization is performed in the frequency range 2 -18 GHz in a multipath environment and in an anechoic chamber, respectively. Fig. 2 shows the measuring setup in a non-anechoic environment (Fig. 2(a)), using two horn antennas with an aperture size of 16.5 by 16.5 cm², 18 and the setup used for measuring inside an anechoic chamber (Fig. 2(b)) with the same broadband horn antennas. The measurements in the anechoic environment have been performed only at one distance where the material is placed in the far field (Fraunhofer) region from the antenna which corresponds to 75 cm at the lowest frequency, 2 GHz.

In a multipath environment, an absorbing material sample can be placed between the antennas in order to reduce the mutual coupling between them. The mutual coupling impacts on the accuracy of the measurements, resulting in an incorrect extraction of the phase. Another method to reduce the coupling between the antennas, [START_REF] Schultz | Focused beam methods: Measuring microwave materials in free space[END_REF] consists in placing the antennas as far as possible one from the other [START_REF] Schultz | Focused beam methods: Measuring microwave materials in free space[END_REF] . The optimal distance between the antennas is 2D 2 /λ at the lowest frequency. This represents the beginning of the far field region. In order to remove any unwanted residual reflections and refractions left after the "deembedding" process, six different distances Antennas -Material between 75 cm (2D 2 /λ at 2 GHz is 72.6 cm) and 125 cm were used [START_REF] Tamas | [END_REF] . At each distance, four transmission coefficients were measured i.e., for the antennas alone, for a metallic plate, for a material sample, and for a material sample with a metallic plate placed on its back. For each distance we apply the "de-embedding" process 7,[START_REF] Hanson | Instrumentation and Measurement[END_REF][START_REF] Balanis | Antenna Theory Analysis and Design Third Edition[END_REF] and we obtain a corresponding reflection coefficient . The average reflection coefficient is then deduced by applying the mean value of all the like this: .

In order to enhance the elimination of the scattering effects of the environment, another technique is applied. In this case, we have moved the measurement configuration to the left and to the right by keeping the same distance between the antennas and the MUT for each distance as shown in Fig. 3. This type of measurement is actually similar to those used for RCS (Radar Cross -Section) measurements and allows reducing the environment's scattering by averaging the mobile echoes 24 .

It is necessary to move the entire measurement setup at least in two positions in order to reduce the effect of the multipath transmissions.

FIG. 3. Reduction of the effect of moving items in a non-anechoic environment: the entire setup is moved but keeping the same distance between the antennas and the sample.

III. RESULTS AND DISCUSSION

A. Teflon

In order to validate the characterization method in a multipath environment, a known dielectric material with very low losses (Teflon 25 sample) was measured.

Based on the radiation pattern and on the first measurement distance, the optimal thickness, length and width of the sample can be found, in order to reduce the effect of the scattering on its edges for broadband studies.

The optimal thickness of the material sample should be at least λ/4 at the lowest frequency i.e., at 2 GHz.

The distance between the antennas and the sample in far field at 2 GHz is 75 cm. The antennas 18 have a -3dB beamwidth of 59 0 (1.03 radians) at 2 GHz, of at least 11.5 0 (0.2 radians) at 10 GHz and of 8 0 (0.14 radians) at 18 GHz, respectively.

The size of the sample can be approximated to the arc length on the wave front. At 10 GHz we have an arc length of 15 cm for the first distance of 75 cm. Since at 2 GHz the size of the sample would have been very big (77.25 cm) and because our technique is used to reduce the effects of the multipath transmission in an ordinary room, we proposed 15 cm for the length and width of the MUT, which is the arc length for the middle frequency, 10 GHz.

The final size of the sample is 15cm x 15cm x 6cm (Large x width x height). With this dimensions we proved that the accurate reflection coefficient can be extracted by applying our mean value method on the measurements at different distances.

Theoretically, the Material Under Test, here Teflon, has a permittivity of 2.1 and a loss tangent of 0.0003 at 10 GHz.

Two types of measurements have been done. One using the MUT alone showing that we have transmission through it and another measurement for the MUT backed by a metallic plate. This second measurement has been performed because it gives us the absorption of the material and also because it is necessary for the permittivity extraction method.

Fig. 4(a) shows the average reflection coefficient obtained from all 30 measurements (6 different antenna-material distances and 5 different displacements of the measurement setup) after applying the "de-embedding" method for the Teflon sample.

Fig. 4(b) shows the reflection coefficient, simulated using 3D CST Microwave Studio and measured for the Teflon sample. We can see that the measured result (dotted red curve) is very different from the one obtained by simulation; when applying the "de-embedding" procedure and the mean value method we can extract an exact result (solid red curve) similar to the simulation and to the measurement in the anechoic chamber.

The phase measured in the multipath environment without "de-embedding" is represented from 10 GHz and is fluctuating dramatically (Fig. 4(c)) over the entire frequency range and therefore it cannot be compared with the theoretical phase. The only reason to depict the phase variation in Fig. 4(c) was to better understand why the "de-embedding" procedure and the mean value are so important for the determination of the correct reflection coefficient of a material. The phase extracted from the same measurement, after the "de-embedding" and applying the average value is confirmed correct by comparison with the measurement in the anechoic chamber as well as from simulation. (see Fig. 4(c)). Now, we consider the Teflon sample backed by a metallic plate (Fig. 5). Fig. 5(a) compares the average values of the reflection coefficient obtained in the real environment with and without "deembedding", with the results obtained in an anechoic chamber and with the CST simulation results.

The phase measured without "de-embedding" and without averaging has the same variation as the one shown in Fig. 4(c). The phase extracted by applying the "de-embedding" and by averaging the

B. Absorbing material

We have also measured a commercial absorbing material named ECCOSORB LS 22 26 . This material is a polyurethane foam loaded with carbon powder. From the material characteristics 26 , we can note that for the frequency range 2 GHz to 18 GHz, the permittivity is varying from 2.6 to 1.4 and the loss tangent is decreasing from 1.8 to 0.6. A sample with a size of 15 cm x 15 cm x 5 cm was characterized.

As before, we measured at several distances four transmission coefficients: with the antennas alone, with a metallic plate, with a material sample, and with a material sample backed with a metallic plate. We extracted from the measured data all the reflection coefficients and the real part of the permittivity for our material. The results from the measurement are close to those achieved from simulation. By measuring the reflection coefficient of the material sample backed by the metallic plate, (Fig. 6(b)) we can find the absorption coefficient and assess whether or not the material can be considered lossy. Nevertheless, we could have an idea on how high the losses are by using an optimization algorithm under CST for the reflection coefficient obtained by measurement. We applied the technique described by Fenner et al. [START_REF] Fenner | [END_REF] to extract the real part of the permittivity from measurement. For the simulation we use this real part of the permittivity and proposed an initial value for the imaginary part of the permittivity that is optimized in CST Microwave Studio. Thus, the optimization finds an imaginary part of the permittivity that gives us the complex ε r used to determinate the reflection coefficient simulated that is as close as possible with the measurement results. The method used for extracting the permittivity is that described by Fenner et al. [START_REF] Fenner | [END_REF]27 . The equations used are shown below:

-TE polarization (5)

where Z 2 is the transverse impedance and Z A and Z B are the wave impedances of the material sample, respectively backed by the metallic plate (MUT -MP) and alone (MUT). Since the measurements have been done in free space, the wave impedance is Z 0 = 377 Ω.

Our materials used are non-magnetic, thus the permeability is μ = μ 0 . This method uses the impedance for the extraction of the material's properties. The phase reference for this technique has to be at the front surface of the material, this is why the "de-embedding" process is very important.

Because we have used a bistatic technique the angle of incidence Ф is 6.31 0 for a distance of 75 cm between the antennas and the material and gets smaller the further we go in far field. Fig. 7 shows the real part of permittivity for a Teflon sample and for the absorbing material, extracted from the measurements in the frequency range 2-18 GHz. For the Teflon sample (Fig. 7(a))

the mean value for the permittivity in this frequency range is 2.06. This result is comparable to the theoretical value (obtained by applying ( 5) -( 8) on the reflection coefficients obtained from simulation) and to the result obtained from the measurement in the anechoic chamber. Note that in CST we used a permittivity of ɛ' = 2.1 at 10 GHz. In CST we are considering an ideal environment with no outside interference and simulate using the sample block as an infinite structure having the same thickness as in the measurement.

Moreover, in Fig. 7(a) we have presented measured results in the X band for a Teflon sample in a waveguide. This result obtained using the Nicolson -Ross -Weir method 28,[START_REF] Weir | Proceedings of the IEEE[END_REF] is close to that obtained in a multipath environment, to the theoretical value in CST, coaxial probe measurement 30 and anechoic chamber. the value given by the fabricant. The results obtained in a waveguide, coaxial probe and anechoic chamber confirm also the exactitude of our method in a real environment.

CONCLUSION

We have developed a method for material characterization in a multipath environment. The proposed measured method does not require an anechoic chamber and allows the extraction of the material's complex reflection coefficient and real part of the permittivity. Thus, measurements can be performed into an ordinary room, by using the distance averaging method, in order to reduce the effects of the multipath propagation. The reflection coefficient that we found for a Teflon sample or a lossy material commercial sample is close to the simulated one, as well as to the reflection coefficient measured in an anechoic chamber. By using this parameter, the real part of the permittivity was obtained. This value is close to the theoretical value and to that obtained in waveguide and anechoic chamber.

FIG. 1 .

 1 FIG. 1. Measurement configuration with the horn antennas.

  FIG. 2. Measurement configurations in (a) a multipath environment and (b) in anechoic chamber.

FIG. 4 .

 4 FIG. 4. Measurement and simulation results of the reflection coefficient (a) mean value in red and the superposition of all measurements (b) magnitude and (c) phase for a Teflon sample.

5 .

 5 chamber Measurement in real environment with "de-embedding" Measurement in real environment without "de-embedding" a multipath environment is close to that obtained in the anechoic chamber and to the simulation as can be seen in Fig.5(b). Measurement and simulation results of the reflection coefficient (a) magnitude and (b) phase for a Teflon sample backed with a metallic plate.

6 .

 6 chamber Measurement in real environment with "de-embedding" Measurement in real environment without "de-embedding" 6(a) we compared the simulation results with the measurement results in the anechoic chamber and in the non-anechoic environment for the absorbing material alone. Measured and simulated results for the reflection coefficient of (a) ECCOSORB LS 22 without a metallic plate and (b) ECCOSORB LS 22 backed with a metallic plate.

  environment w ith "de-embedding" Measurement in real environment w ithout "de-embedding" Simulation CST using the material characteristic given in [26] environment w ith "de-embedding" Measurement in real environment w ithout "de-embedding" Simulation CST using the material characteristic given in[26]

FIG. 7 .

 7 FIG. 7. Permittivity of (a) Teflon and (b) ECCOSORB LS 22 loaded with carbon powder.Fig. 7(b) presents the measured results for the real part of the permittivity of the loaded

ACKNOWLEDGMENTS

The authors would like to thank the Direction Générale de l'Armement (DGA) for its financial support.