Thomas Carle

Manel Djemal

Dumitru Potop-Butucaru

Robert De Simone

Zhen Zhang

Static mapping of real-time applications onto massively parallel processor arrays

On-chip networks (NoCs) used in multiprocessor systems-on-chips (MPSoCs) pose significant challenges to both on-line (dynamic) and off-line (static) real-time scheduling approaches. They have large numbers of potential contention points, have limited internal buffering capabilities, and network control operates at the scale of small data packets. Therefore, efficient resource allocation requires requires scalable algorithms working on hardware models with a level of detail that is unprecedented in real-time scheduling. We consider here a static scheduling approach, and we target massively parallel processor arrays (MPPAs), which are MPSoCs with large numbers (hundreds) of processing cores. We first identify and compare the hardware mechanisms supporting precise timing analysis and efficient resource allocation in existing MPPA platforms. We determine that the NoC should ideally provide the means of enforcing a global communications schedule that is computed off-line (before execution) and which is synchronized with the scheduling of computations on processors. On the software side, we propose a novel allocation and scheduling method capable of synthesizing such global computation and communication schedules covering all the execution, communication, and memory resources in an MPPA. To allow an efficient use of the hardware resources, our method takes into account the specificities of MPPA hardware and implements advanced scheduling techniques such as software pipelining and pre-computed preemption of data transmissions. We evaluate our technique by mapping two signal processing applications, for which we obtain good latency, throughput, and resource use figures.

I. INTRODUCTION

One crucial problem in real-time scheduling is that of ensuring that the application software and the implementation platform satisfy the hypotheses allowing the application of specific schedulability analysis techniques [START_REF] Layland | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF]. Such hypotheses are the availability of a priority-driven scheduler or the possibility of including scheduler-related costs in the durations of the tasks. Classical work on real-time scheduling [START_REF] Davis | A survey of hard real-time scheduling for multiprocessor systems[END_REF], [START_REF] Goossens | Priority-driven scheduling of periodic task systems on multiprocessors[END_REF], [START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF] has proposed formal models allowing schedulability analysis in classical mono-processor and distributed settings. But the advent of multiprocessor systems-on-chip (MPSoC) architectures imposes significant changes to these models and to the scheduling techniques using them.

MPSoCs are becoming prevalent in both general-purpose and embedded systems. Their adoption is driven by scalable performance arguments (concerning speed, power, etc.), but this scalability comes at the price of increased complexity of both the software and the software mapping (allocation and scheduling) process.

Part of this complexity can be attributed to the steady increase in the quantity of software that is run by a single system. But there are also significant qualitative changes concerning both the software and the hardware. In software, more and more applications include parallel versions of classical signal or image processing algorithms [START_REF] Villalpando | Investigation of the tilera processor for real time hazard detection and avoidance on the altair lunar lander[END_REF], [START_REF] Aubry | Extended cyclostatic dataflow program compilation and execution for an integrated manycore processor[END_REF], [21], which are best modeled using data-flow models (as opposed to independent tasks). Providing functional and real-time correctness guarantees for parallel code requires an accurate control of the interferences due to concurrent use of communication resources. Depending on the hardware and software architecture, this can be very difficult [START_REF] Wilhelm | Embedded systems: Many cores -many problems (invited paper)[END_REF], [START_REF] Herlihy | The Art of Multiprocessor Programmin[END_REF].

Significant changes also concern the execution platforms, where the gains predicted by Moore's law no longer translate into improved single-processor performance, but in a rapid increase of the number of processor cores placed on a single chip [START_REF] Borkar | Thousand core chips -a technology perspective[END_REF]. This trend is best illustrated by the massively parallel processor arrays (MPPAs), which we target in this paper. MPPAs are MPSoCs characterized by:

• Large numbers of processing cores, ranging in current silicon implementations to a few tens to a few hundreds [START_REF]The TilePro64 many-core architecture[END_REF], [START_REF]The MPPA256 many-core architecture[END_REF], [START_REF]The Epiphany many-core architecture[END_REF], [START_REF]A 48-core ia-32 processor in 45nm cmos using on-die message-passing and dvfs for performance and power scaling[END_REF]. The cores are typically chosen for their area or energy efficiency instead of raw computing power.

• A regular internal structure where processor cores and internal storage (RAM banks) are divided among a set of identical tiles, which are connected through one or more NoCs with regular structure (e.g. torus, mesh).

Industrial [START_REF]The TilePro64 many-core architecture[END_REF], [START_REF]The MPPA256 many-core architecture[END_REF], [START_REF]The Epiphany many-core architecture[END_REF], [START_REF]A 48-core ia-32 processor in 45nm cmos using on-die message-passing and dvfs for performance and power scaling[END_REF] and academic [START_REF] Goossens | AEthereal network on chip: Concepts, architectures, and implementations[END_REF], [START_REF] Bjerregaard | Implementation of guaranteed services in the mango clockless network-on-chip[END_REF], [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF] MPPA architectures targeting hard real-time applications already exist, but the problem of mapping applications on them remains largely open. There are two main reasons to this. The first one concerns the NoCs: as the tasks are more tightly coupled and the number of resources in the system increases, the on-chip networks become critical resources, which need to be explicitly considered and managed during real-time scheduling. Recent work [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF], [START_REF] Kashif | Ortap: An offset-based response time analysis for a pipelined communication resource model[END_REF], [START_REF] Nikolic | Are virtual channels the bottleneck of priority-aware wormhole-switched noc-based manycores?[END_REF] has determined that NoCs have distinctive traits requiring significant changes to classical real-time scheduling theory [START_REF] Goossens | Priority-driven scheduling of periodic task systems on multiprocessors[END_REF]. The second reason concerns automation: the complexity of MPPAs and of the (parallel) applications mapped on them is such that the allocation and scheduling must be largely automated.

We address here this automation need, and unlike previous work on the subject we focus on static (off-line) scheduling approaches. In theory, off-line algorithms allow the computation of scheduling tables specifying an optimal allocation and real-time scheduling of the various computations and communications onto the resources of the MPPA. In practice, this ability is severely limited by 3 factors:

1) The application may exhibit a high degree of dynamicity due to either environment constraints or to execution time variability resulting from data-dependent conditional control. 12) The hardware may not allow the implementation of optimal scheduling tables. For instance, most MPPA architectures provide only limited control over the scheduling of communications inside the NoC.

3) The mapping problems we consider are NP-hard. In practice, this means that optimality cannot be attained, and that efficient heuristics are needed. In the remainder of the paper we assume that our applications are static enough to benefit from off-line scheduling (e.g. parallelized versions of signal processing and embedded control algorithms).

Our contributions concern the hardware and the mapping. We start with an in-depth review of MPPA/NoC architectures with support for real-time scheduling. This analysis allows us to determine that NoCs allowing static communication scheduling offer the best support to off-line application mapping. This analysis also allows us to define a new technique and tool, called LoPhT, for automatic real-time mapping and code generation, whose global flow is pictured in Fig. 1. and precise hardware descriptions including all potential NoC contention points. It uses advanced off-line scheduling techniques such as software pipelining and pre-computed preemption, and it takes into account the specificities of the MPPA hardware to build scheduling tables that provide good latency and throughput guarantees and ensure an efficient use of computation and communication resources. Scheduling tables are then automatically converted into sequential code ensuring the correct ordering of operations on each resource and the respect of the real-time guarantees. We provide results for two applications that are mapped onto the MPPA platform described in [START_REF] Djemal | Programmable routers for efficient mapping of applications onto NoC-based MPSoCs[END_REF]. Our results show that the off-line mapping of communications not only allows us to provide static latency and throughput guarantees, but may also improve the speed of the application.

Our tool takes as input data-flow synchronous specifications

II. MPPA/NOC ARCHITECTURES FOR THE REAL-TIME

We start with a general introduction to MPPA platforms, and then present the main characteristics of existing NoCs. We are mainly interested here in the traffic management mechanisms supporting real-time implementation.

A. Structure of an MPPA

Our work concerns hard real-time systems where timing guarantees must be determined by static analysis methods before system execution. Complex memory hierarchies involving multiple cache levels and cache coherency mechanisms are known to complicate timing analysis [START_REF] Wilhelm | Embedded systems: Many cores -many problems (invited paper)[END_REF], [START_REF] Hardy | Wcet analysis of multi-level non-inclusive setassociative instruction caches[END_REF], and we assume they are not used in the MPPA platforms we consider. Under this hypothesis, all data transfers between tiles are performed through one or more NoCs.

A NoC can be described in terms of point-to-point communication links and NoC routers which perform the routing and scheduling (arbitration) functions. Fig. 2 provides the description of a 2-dimensional (2D) mesh NoC like the ones used in the Adapteva Epiphany [START_REF]The Epiphany many-core architecture[END_REF], Tilera TilePro [START_REF]The TilePro64 many-core architecture[END_REF], or DSPIN [START_REF] Panades | A low cost networkon-chip with guaranteed service well suited to the GALS approach[END_REF]. The structure of a router in a 2D mesh NoC is 3), and a scheduling/arbitration component, which we call multiplexer (labeled M). Data enters the router through demultiplexers and exits through multiplexers.

To allow on-chip implementation at a reasonable area cost, NoCs use very simple routing algorithms. For instance, the Adapteva [START_REF]The Epiphany many-core architecture[END_REF], Tilera [START_REF]The TilePro64 many-core architecture[END_REF], and DSPIN NoCs [START_REF] Panades | A low cost networkon-chip with guaranteed service well suited to the GALS approach[END_REF] use an Xfirst routing algorithm where data first travels all the way in the X direction, and only then in the Y direction. Furthermore, all NoCs mentioned in this paper use simple wormhole switching approaches [START_REF] Ni | A survey of wormhole routing techniques in direct networks[END_REF] requiring that all data of a communication unit (e.g. packet) follow the same route, in order, and that the communication unit is not logically split during transmission. The use of a wormhole switching approach is justified by the limited buffering capabilities of NoCs [START_REF] Nikolic | Are virtual channels the bottleneck of priority-aware wormhole-switched noc-based manycores?[END_REF] and by the possibility of decreasing transmission latencies (by comparison with more classical store-and-forward approaches). But the use of wormhole switching means that one data transmission unit (such as a packet) is seldom stored in a single router buffer. Instead, a packet usually spans over several routers, so that its transmission strongly synchronizes multiplexers and demultiplexers along its path.

B. Support for real-time implementation

Given the large number of potential contention points (router multiplexers), and the synchronizations induced by data transmissions, providing tight static timing guarantees is only possible if some form of flow control mechanism is used.

In NoCs based on circuit switching [START_REF]Pnoc: a flexible circuit-switched noc for fpga-based systems[END_REF], inter-tile communications are performed through dedicated communication channels formed of point-to-point physical links. Two channels cannot share a physical link. This is achieved by statically fixing the output direction of each demultiplexer and the data source of each multiplexer along the channel path. Timing interferences between channels are impossible, which radically simplifies timing analysis, and the latency of communications is low. But the absence of resource sharing is also the main drawback of circuit switching, resulting in low numbers of possible communication channels and low utilization of the NoC resources. Reconfiguration is usually possible, but it carries a large timing penalty.

Virtual circuit switching is an evolution of circuit switching which allows resource sharing between circuits. But resource sharing implies the need for arbitration mechanisms inside NoC multiplexers. Very interesting from the point of view of timing predictability are NoCs where arbitration is based on time division multiplexing (TDM), such as Aethereal [START_REF] Goossens | AEthereal network on chip: Concepts, architectures, and implementations[END_REF], Nostrum [START_REF] Millberg | Guaranteed bandwidth using looped containers in temporally disjoint networks within the Nostrum network-on-chip[END_REF], and others [START_REF] Sorensen | A light-weight statically scheduled network-on-chip[END_REF]. In a TDM NoC, all routers share a common time base. The point-to-point links are reserved for the use of the virtual circuits following a fixed cyclic schedule (a scheduling table). The reservations made on the various links ensure that communications can follow their path without waiting. TDM-based NoCs allow the computation of precise latency and throughput guarantees. They also ensure a strong temporal isolation between virtual circuits, so that changes to a virtual circuit do not modify the real-time characteristics of the other.

When no global time base exists, the same type of latency and throughput guarantees can be obtained in NoCs relying on bandwidth management mechanisms such as Kalray MPPA [START_REF]The MPPA256 many-core architecture[END_REF], [START_REF] Harrand | Network on chip with quality of service[END_REF]. The idea here is to ensure that the throughput of each virtual circuit is limited to a fraction of the transmission capacity of a physical point-to-point link, by either the emitting tile or by the NoC routers. Two or more virtual circuits can share a point-to-point link if their combined transmission needs are less than what the physical link provides.

But TDM and bandwidth management NoCs have certain limitations: One of them is that latency and throughput are correlated [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF], which may result in a waste of resources. But the latency-throughput correlation is just one consequence of a more profound limitation: TDM and bandwith management NoCs largely ignore the fact that the needs of an application may change during execution, depending on its state. For instance, when scheduling a data-flow synchronous graph with the objective of reducing the duration of one computation cycle (also known as makespan or latency), it is often useful to allow some communications to use 100% of the physical link, so that they complete faster, before allowing all other communications to be performed.

One way of taking into account the application state is by using NoCs with support for priority-based scheduling [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF], [START_REF] Panades | A low cost networkon-chip with guaranteed service well suited to the GALS approach[END_REF], [START_REF]A 48-core ia-32 processor in 45nm cmos using on-die message-passing and dvfs for performance and power scaling[END_REF]. In these NoCs, each data packet is assigned a priority level (a small integer), and NoC routers allow higherpriority packets to pass before lower-priority packets. To avoid priority inversion phenomenons, higher-priority packets have the right to preempt the transmission of lower-priority ones. In turn, this requires the use of one separate buffer for each priority level in each router multiplexer, a mechanism known as virtual channels (VCs) in the NoC community [START_REF] Panades | A low cost networkon-chip with guaranteed service well suited to the GALS approach[END_REF].

The need for VCs is the main limiting factor of prioritybased arbitration in NoCs. Indeed, adding a VC is as complex as adding a whole new NoC [START_REF] Yoon | Virtual channels vs. multiple physical networks: a comparative analysis[END_REF], [START_REF] Carara | Router architecture for highperformance nocs[END_REF], and NoC resources (especially buffers) are expensive in both power consumption and area [START_REF] Moscibroda | A case for bufferless routing in on-chip networks[END_REF]. To our best knowledge, among existing silicon implementations only the Intel SCC chip offers a relatively large numbers of VCs (eight) [START_REF]A 48-core ia-32 processor in 45nm cmos using on-die message-passing and dvfs for performance and power scaling[END_REF], and it is targeted at high-performance computing applications. Industrial MPPA chips targeting an embedded market usually feature multiple, specialized NoCs [START_REF]The TilePro64 many-core architecture[END_REF], [START_REF]The Epiphany many-core architecture[END_REF], [START_REF] Harrand | Network on chip with quality of service[END_REF] without virtual channels. Other NoC architectures feature low numbers of VCs. Current research on priority-based communication scheduling has already integrated this limitation, by investigating the sharing of priority levels [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF].

Significant work already exists on the mapping of real-time applications onto priority-based NoCs [START_REF] Shi | Schedulability analysis and task mapping for real-time on-chip communication[END_REF], [START_REF] Racu | Using genetic algorithms to map hard realtime on noc-based systems[END_REF], [START_REF] Nikolic | Are virtual channels the bottleneck of priority-aware wormhole-switched noc-based manycores?[END_REF], [START_REF] Kashif | Ortap: An offset-based response time analysis for a pipelined communication resource model[END_REF]. This work has shown that priority-based NoCs support the efficient mapping of independent tasks.

But we already explained that the large number of computing cores in an MPPA means that applications are also likely to include parallelized code which is best modeled by large sets of relatively small dependent tasks (data-flow graphs) with predictable functional and temporal behavior [START_REF] Villalpando | Investigation of the tilera processor for real time hazard detection and avoidance on the altair lunar lander[END_REF], [START_REF] Aubry | Extended cyclostatic dataflow program compilation and execution for an integrated manycore processor[END_REF], [21]. Such timing-predictable specifications are those that can a priori take advantage of a static scheduling approach, which provides best results on architectures with support for static communication scheduling [START_REF]The TilePro64 many-core architecture[END_REF], [START_REF]Baring it all to software: The raw machine[END_REF], [START_REF] Djemal | Programmable routers for efficient mapping of applications onto NoC-based MPSoCs[END_REF]. Such architectures allow the construction of an efficient (possibly optimal) global computation and communication schedule, represented with a scheduling table and implemented as a set of synchronized sequential computation and communication programs. Computation programs run on processor cores to sequence task executions and the initiation of communications. Communication programs run on specially-designed microcontrollers that control each NoC multiplexer to fix the order in which individual data packets are transmitted. Synchronization between the programs is ensured by the data packet communications themselves.

Like in TDM NoCs, the use of global computation and communication scheduling tables allows the computation of very precise latency and throughput estimations. Unlike in TDM NoCs, NoC resource reservations can depend on the application state. Global time synchronization is not needed, and existing NoCs based on static communication scheduling do not use it [START_REF]The TilePro64 many-core architecture[END_REF], [START_REF]Baring it all to software: The raw machine[END_REF], [START_REF] Djemal | Programmable routers for efficient mapping of applications onto NoC-based MPSoCs[END_REF]. Instead, global synchronization is realized by the data transmissions (which eliminates some of the run-time pessimism of TDM-based approaches).

The microcontrollers that drive each NoC router multiplexer are similar in structure to those used in TDM NoCs to enforce the TDM reservation pattern. The main difference is that the communication programs are usually longer than the TDM configurations, because they must cover longer execution patterns. This requires the use of larger program memory (which can be seen as part of the tile program memory [START_REF] Djemal | Programmable routers for efficient mapping of applications onto NoC-based MPSoCs[END_REF]). But like in TDM NoCs, buffering needs are limited and no virtual channel mechanism is needed, which results in lower power consumption.

From a mapping-oriented point of view, determining exact packet transmission orders cannot be separated from the larger problem of building a global scheduling table comprising both computations and communications. By comparison, mapping onto MPPAs with TDM-based or bandwith reservation-based NoCs usually separates task allocation and scheduling from the synthesis of a NoC configuration independent from the application state [START_REF] Lu | Tdm virtual-circuit configuration for networkon-chip[END_REF], [START_REF] Aubry | Extended cyclostatic dataflow program compilation and execution for an integrated manycore processor[END_REF], [START_REF] Zhai | Exploiting just-enough parallelism when mapping streaming applications in hard real-time systems[END_REF] .

Under static communication scheduling, there is little runtime flexibility, as all scheduling possibilities must be considered during the off-line construction of the global scheduling table. For very dynamic applications this can be difficult. This is why existing MPPA architectures that allow static communication scheduling also allow communications with dynamic (Round Robin) arbitration.

In conclusion, NoCs allowing static communication scheduling offer the best temporal precision in the off-line mapping of dependent tasks (data-flow graphs), while prioritybased NoCs are better at dealing with more dynamic applications. As future systems will include both statically parallelized code and more dynamic aspects, NoCs will include mechanisms supporting both off-line and on-line communication scheduling. Significant work already exists on the real-time mapping for priority-based platforms, while little work has addressed the NoCs with static communication scheduling. In the remainder of the paper we address this issue by proposing a new off-line scheduling technique.

III. RELATED WORK

The previous section already discussed existing work on taking into account NoC structure during real-time mapping.

But our work is also closely related to classical results on the off-line real-time mapping of dependent task systems onto multiprocessor and distributed architectures [START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF], [START_REF] Xu | Multiprocessor scheduling of processes with release times, deadlines, precedence, and exclusion relations[END_REF], [START_REF] Eles | Scheduling with bus access optimization for distributed embedded systems[END_REF], [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives: a seamless flow of graphs transformations[END_REF]. The objective is the same as in our case: the synthesis of optimized time-triggered execution patterns. But there are also significant differences. Most NoCs rely on wormhole routing, which requires synchronized reservation of the resources along the communication paths. By comparison, the cited papers either use a store-and-forward routing paradigm [START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF], [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives: a seamless flow of graphs transformations[END_REF], [START_REF] Eles | Scheduling with bus access optimization for distributed embedded systems[END_REF] that is inapplicable to NoCs, or simply do not model the communication media [START_REF] Xu | Multiprocessor scheduling of processes with release times, deadlines, precedence, and exclusion relations[END_REF].

A second difference concerns the complexity of the architecture description. MPPAs have more computation and communication resources than typical distributed architectures considered in classical real-time. Moreover, resource characterization has clock cycle precision. To scale up without losing timing precision, we employ scheduling heuristics of low computational complexity, avoiding the use of techniques such as backtracking [START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF], but taking advantage of low-level architectural detail concerning the NoC.

A third difference concerns fine-grain resource allocation and scheduling operations, such as the allocation of data variables into memory banks or the scheduling of DMA commands onto the processor cores. These operations are often overlooked in classical distributed scheduling. But on an MPPA platform they can have a major impact on execution durations. This why we need to explicitly consider them at off-line mapping time.

More generally, our work is also close to previous results on static, but not real-time, mapping of applications onto MPPAs. The StreamIt compiler [START_REF] Amarasinghe | Language and compiler design for streaming applications[END_REF] uses scheduling tables as an internal compiler representation. But StreamIt uses timing information coming from simulations to guide the mapping process, and does not take into account timing interferences due to the mapping itself. Previous work on the compilation of the SigmaC language onto the Kalray MPPA platform [START_REF] Aubry | Extended cyclostatic dataflow program compilation and execution for an integrated manycore processor[END_REF] also uses timing information coming from simulations, but its objective is to generate code for a micro-kernel employing dynamic task scheduling.

Similar approaches are taken in more dynamic techniques aimed at signal processing systems [START_REF] Bhattacharyya | Handbook of Signal Processing Systems[END_REF], [START_REF] Genius | Space optimal solution for data reordering in streaming applications on noc based mpsoc[END_REF], [START_REF] Bebelis | A framework to schedule parametric dataflow applications on many-core platforms[END_REF]. Techniques such as DOL [START_REF] Bacivarov | Methods and tools for mapping process networks onto multi-processor systems-on-chip[END_REF] or the one of Zhai et al. [START_REF] Zhai | Exploiting just-enough parallelism when mapping streaming applications in hard real-time systems[END_REF] allow for real-time mapping, but without considering the details of the NoC. There is a price to pay for this: In DOL, mapping results are not proved safe, because timing information comes from simulations which may not cover the worst case. In the approach of Zhai et al., NoC resource allocation that guarantees the latency of inter-tile communications is done before application mapping, and is not optimized afterwards. Our method of taking into account low-level NoC details could potentially be used to improve such mapping techniques.

IV. HARDWARE MODEL

A. MPPA hardware

We use in this paper the MPPA architecture of Djemal et al. [START_REF] Djemal | Programmable routers for efficient mapping of applications onto NoC-based MPSoCs[END_REF]. The overall architecture of the MPPA and of its routers is that of Figures 2 and3. NoC control operates at packet level, each packet being a bounded sequence of flits (FLow control unITs, the data unit that can be sent over a NoC link in one clock cycle). NoC routing if of X-first type.

Packet scheduling (arbitration) in the router multiplexers can be either dynamic (with a Round Robin policy), or fixed off-line. In the second case, the order in which packets are transmitted is specified with a communication program stored in the router controller component of Fig. 4. Regardless of the arbitration type, the transmission of a packet can never be interrupted. The NoC only allows the buffering of 3 flits for each link. For efficiency and timing predictability purposes, the tiles of our MPPA are formed of the following components (Fig. 4):

• 16 MIPS32 processor cores with separate instruction and data caches.

• 32 data memory banks that can be accesses independently, for a total of maximum 2 Mbytes of RAM. The program is stored in a separate RAM or ROM bank.

• A DMA unit which allows the queuing of multiple DMA commands in order to diminish timing interference between computations and communications.

• A hardware lock component. The locks allow lowoverhead synchronization and make it easier to preserve temporal predictability because they do not alter the cache state (as opposed to interrupt-based synchronization).

Locks are used to ensure exclusive access to data RAM banks.

• A NoC router controller that runs the arbitration programs of the NoC router.

All these components are linked together through a full crossbar local interconnect which is also linked to the NoC.

B. Resource modeling

To allow off-line mapping onto our architectures, we need to identify the set of abstract computation and communication resources that are considered during allocation and scheduling.

We associate one communication resource to each of the multiplexers of NoC routers and to each DMA. We name them as follows: N (i, j)(k, l) is the inter-router wire going from T ile(i, j) to T ile(k, l); In(i, j) is the output of the router (i, j) to the local tile; DM A(i, j) is the output of T ile(i, j) to the local router.

This paper focuses on NoC modeling and handling. To this end, we consider in this paper a resource model that simplifies as much as possible the representation of the computing tiles. All the 16 processor cores of the tile are seen as a single, very fast computing resource. This means that operations will be allocated to the tile as if it were a sequential processor, but the allocated operations are in fact parallel code running on all 16 processors. In Fig. 2 there are just 12 tile resources representing 192 processor cores. This simplification largely reduces the complexity of our presentation, and also satisfies our evaluation needs, given that the 2 applications used as examples can be organized into operations that are easily parallelized onto 16 processors.

C. Communication durations

All inter-tile data transmissions are performed using the DMA units. If a transmission is not blocked on the NoC, then its duration on the sender side only depends on the size of the transmitted data. The exact formula is d = s + ⌈s/maxpayload⌉ * PacketHeaderSize, where d is the duration in clock cycles of the DMA transfer from the start of the transmission to the cycle where a new transmission can start, s is the data size in 32-bit words, MaxPayload is the maximum payload of a NoC packet produced by the DMA (in 32-bit words), and PacketHeaderSize is the number of cycles that are lost for each packet in the chosen NoC. In our case, MaxPayload=16 flits and PacketHeaderSize=4 flits.

In addition to this transmission duration, we must also account in our computations for:

• The DMA transfer initiation, which consists in 3 uncached RAM accesses plus the duration of the DMA reading the payload of the first packet from the data RAM. This cost is over-approximated as 30 cycles.

• The latency of the NoC, which is the time needed for one flit to traverse the path from source to destination. This latency is of 3 * n, where n is the number of NoC multiplexers on the route of the transmission.

V. APPLICATION SPECIFICATION

The application to be mapped is specified under the form of a data-flow graph (known as a dependent task system in the real-time community). The input of our tool is a fullfledged data-flow synchronous formalism like that of [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives: a seamless flow of graphs transformations[END_REF]. For presentation reasons, however, we simplify it as follows.

A data-flow graph D is a directed graph with two types of arcs D = {T (D), A(D), ∆(D)}. Here, T (D) is the finite set of tasks (data-flow blocks). The finite set A(D) contains dependencies of the form a = (src(a), dst(a), type(a)), where src(a), dst(a) ∈ T (D) are the source, respectively the destination task of a, and type(a) is the type of data transmitted from src(a) to dst(a). The directed graph determined by A(D) must be acyclic. The finite set ∆(D) contains delayed dependencies of the form δ = (src(δ), dst(δ), type(δ), depth(δ)), where src(δ), dst(δ), type(δ) have the same meaning as for simple dependencies and depth(δ) is a strictly positive integer called the depth of the dependency.

Data-flow graphs have a cyclic execution model. At each execution cycle of the task set, each of the tasks is executed exactly once. We denote with t n the instance of task t ∈ T (D) for cycle n. The execution of the tasks inside a cycle is partially ordered by the dependencies of A(D). If a ∈ A(D) then the execution of src(a) n must be finished before the start of dst(a)

n , for all n.

The dependencies of ∆(D) impose an order between tasks of successive execution cycles. If δ ∈ ∆(D) then the execution of src(δ) n must complete before the start of dst(δ) n+depth(δ) , for all n. We exemplify our formalism with the data-flow graph of Fig. 5, which is a simplified version of an automotive platooning application [START_REF] Pradalier | The CyCab: a car-like robot navigating autonomously and safely among pedestrians[END_REF]. In our figure, each block is a task, solid arcs are simple dependencies, and the dashed arc is a delayed dependency of depth 2. The application is run by a car to determine the position (distance and angle) of another car moving in front of it. It works by cyclically capturing an input image of fixed size. This image is passed through an edge-detecting Sobel filter and then through a histogram search to detect dominant edges (horizontal and vertical). This information is used by the detection and correction function to determine the position of the front car. The whole process is monitored on a display. The delayed dependency represents a feedback from the detection and correction function that allows the adjustment of image capture parameters.

The Sobel filter and the histogram search are parallelized. Each of the Sobel H and Sobel V functions receives one sixth of the whole image (a horizontal slice).

A. Non-functional specification

For each task τ ∈ T (D), we define W CET (τ ) to be a safe upper bound for the worst-case execution time (WCET) of τ on an MPPA tile, in isolation. Note that the WCET values we require are for parallel code running on all the 16 processors of a tile. Tight WCET bounds for such code can be computed using the analysis technique proposed in [START_REF] Puaut | Integrated worst-case execution time estimation of multicore applications[END_REF].

For each data type t associated with a dependency (simple or delayed), we define the worst-case memory footprint of a value of type t. This information allows the computation of the worst-case communication time (WCCT) for a data of that type, using the formula of Section IV-C.

Allocation constraints specify on which tiles a given dataflow block can be executed. In our example, they force the allocation of the capture and display functions onto specific MPPA tiles. More generally, they can be used to confine an application to part of the MPPA, leaving the other tiles free to execute other applications.

VI. OFF-LINE MAPPING AND CODE GENERATION

A. The problem

The real-time mapping and code generation problem we consider in this paper is a bi-criteria optimization problem: Given a data-flow graph and a non-functional specification, synthesize a real-time implementation that minimizes execution cycle latency (duration) and maximizes throughput2 , with priority given to latency. We chose this scheduling problem because it is meaningful in embedded systems design and because its simple definition allows us to focus on the handling of NoC-related issues. Variations of our mapping algorithms can also handle periodicity, start date, and deadline requirements [START_REF] Carle | From dataflow specification to multiprocessor partitioned time-triggered realtime implementation[END_REF].

Our allocation and scheduling problem being NP-complete, we do not aim for optimality. Instead, we rely on lowcomplexity heuristics that allow us to handle large numbers of resources with high temporal precision. Mapping and code generation is realized in 3 steps: The first step produces a latency-optimizing scheduling table using the algorithms of Section VI-B. The second step uses the software pipelining algorithms of [12] (not detailed here) to improve throughput while not changing latency. Finally, once a scheduling table is computed, it is implemented it in a way that preserves its real-time properties, as explained in Section VI-C

B. Latency-optimizing scheduling routine

Our scheduling routine builds a global scheduling table covering all MPPA resources. It uses a non-preemptive scheduling model for the tasks, and a preemptive one for the NoC communications. 3 For each task it reserves exactly one time interval on one of the tiles. For every dependency between two tasks allocated on different tiles, the scheduling routine reserves one or more time intervals on each resource along the route between the two tiles, starting with the DMA of the source tile, and continuing with the NoC multiplexers (recall that the route is fixed under the X-first routing policy).

The scheduling algorithm uses a simple list scheduling heuristic. The tasks of the dependent task system are traversed one by one in an order compatible with the dependencies between them (only the simple dependencies, not the delayed ones, which are considered during throughput optimization). Resource allocation for a task and for all communications associated with its input dependencies is performed upon traversal of the task, and never changed afterwards. Scheduling starts with an empty scheduling table which is incrementally filled as the tasks and the associated communications are reserved time intervals on the various resources.

For each task, scheduling is attempted on all the tiles that can execute the block (as specified by the allocation constraints), at the earliest date possible. Among the possible allocations, we retain the one that minimizes a cost function. This cost function should be chosen so that the final length of the scheduling table is minimized (this length gives the execution cycle latency). Our choice of cost function combines the end date of the task in the schedule (with 95% weight) and the maximum occupation of a CPU in the current scheduling table (with 5% weight). Experience showed that this function produces shorter scheduling tables than the cost function based on task end date alone (as used in [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives: a seamless flow of graphs transformations[END_REF]).

1) Mapping NoC communications:

The most delicate part of our scheduling routine is the communication mapping function MapCommunicationOnPath. When a task is mapped on a tile, this function is called once for each of the input dependencies of the task, if the dependency source is on another tile and if the associated data has not already been transmitted.

Fig. 6 presents a (partial) scheduling table produced by our mapping routine. We shall use this example to give a better intuition on the functioning of our algorithms. We assume here that the execution of task f produces data x which will be used by task g. Our scheduling table shows the result of the mapping of task g onto Tile(2, 2) (which also requires the mapping of the transmission of x) under the assumption that all other tasks (f , h) and data transmissions (y, z, u) were already mapped as pictured (reservations made during the mapping of g have a lighter color). 3 Task preemptions would introduce important temporal imprecision (through the use of interrupts), and are avoided. Data communications over the NoC are naturally divided into packets that are individually scheduled by the NoC multiplexer programs, allowing a form of pre-computed preemption. As part of the mapping of g onto Tile(2, 2), function Map-CommunicationOnPath is called to perform the mapping of the communication of x from Tile(1, 1) to Tile [START_REF] Amarasinghe | Language and compiler design for streaming applications[END_REF][START_REF] Amarasinghe | Language and compiler design for streaming applications[END_REF]. The parameters of its call are Path, StartDate, and DataSize. Parameter Path is the list formed of resources DMA(1, 1), N(1, 1)(1, 2), N(1, 2)(2, 2), and In(2, 2) (the transmission route of x under the X-first routing protocol). Parameter StartDate is given by the end date of task f (in our case 500), and DataSize is the worst-case size of the data associated with the data dependency (in our case 500 32-bit words). Time is measured in clock cycles.

To minimize the overall time reserved for a data transmission, we shall require that it is never blocked waiting for a NoC resource. For instance, if the communication of x starts on the N(1, 1)(1, 2) at date t, then on N(1, 2)(2, 2) it must start at date t + SegmentBufferSize, where SegmentBufferSize is a platform constant defining the time needed for a flit to traverse one NoC resource. In our NoC this constant is 3 clock cycles (in Fig. 6 we use a far larger value of 100 cycles, for clarity).

Building such synchronized rezervation patterns along the communication routes is what function MapCommunica-tionOnPath does. It starts by obtaining the lists of free time intervals of each resource along the communication path, and realigning them by subtracting (i-1) * SegmentBufferSize from the start dates of all the free intervals of the i th resource, for all i. Once this realignment is done on each resource by function ShiftLeftIntervals, finding a reservation along the communication path amounts to finding time intervals that are unused on all resources. To do this, we start by performing (in line 6 of function MapCommunicationOnPath) an intersection operation returning all realigned time intervals that are free on all resources. In Fig. 6, this intersection operation produces (prior to the mapping of x) the intervals [800,1100) and [1400,2100]. The value 2100 corresponds here to the length of the scheduling table prior to the mapping of g.

We then call function ReserveIntervals twice, to make reservations for the data transmission and for the lock command associated with each communication. These two functions produce a list of reserved intervals, which then need to be realigned on each resource. In Fig. 6, these 2 calls reserve the intervals [800,1100), [1400,1700), and [1700,1704). The first 2 intervals are needed for the data transmission, and the third is used for the lock command packet.

2) Multiple reservations: Communications are reserved at the earliest possible date, and function ReserveIntervals allows the fragmentation of a data transmission to allow a better use of NoC resources. In our example, fragmentation allows us to transmit part of x before the reservation for u. If fragmentation were not possible, the transmission of x should be started later, thus delaying the start of g, potentially lengthening the reservation table.

Fragmentation is subject to restrictions arising from the way communications are packetized. More precisely, an interval cannot be reserved unless it has a minimal size, allowing the transmission of at least a packet containing some payload data.

Function ReserveIntervals performs the complex translation from data sizes to needed packets and intervals reservations. We present here an unoptimized version that facilitates understanding. This version reserves one packet at a time, using a free interval as soon as it has the needed minimal size. Packets are reserved until the required DataSize is covered. Like for tasks, reservations are made as early as possible. For each packet reservation the cost of NoC control (under the form of the PacketHeaderSize) must be taken into account.

When the current scheduling table does not allow the mapping of a data communication, function ReserveIntervals will lengthen it so that mapping is possible.

C. Automatic code generation

Once the scheduling table has been computed, executable code is automatically generated as follows: One sequential execution thread is generated for each tile and for each NoC ival := GetFirstInterval(NewIntervalList); 5:

NewIntervalList := RemoveFirstInterval(NewIntervalList); 6:

if IntervalEnd(ival)==ScheduleLength then multiplexer (resources T ile(i, j), N(i, j)(k, l), and In(i, j) in our platform model of Section IV-B). The code of each thread is an infinite loop that executes the (computation or communication) operations scheduled on the associated resource in the order prescribed by their reservations. Recall that each tile contains 16 processor cores, but is reserved as a single sequential resource, parallelism being hidden inside the dataflow blocks. The sequential thread of a tile runs on the first processor core of the tile, but the code of each task can use all 16 processor cores. The code of the NoC multiplexers is executed on the router controllers.

No separate thread is generated for the DMA resource of a tile. Instead, its operations are initiated by thread of the tile. This is possible because the DMA allows the queuing of DMA commands and because mapping is performed so that activation conditions for DMA operations can be computed by the tile resource at the end of data-flow operations. For instance, in the example of Fig. 6, if no other operations are allocated on Tile(0, 0), the two DMA operations needed to send x are queued at the end of f . The synchronization of the threads is realized by explicit lock manipulations by the processors and by the NoC control programs, which force message passing order and implicitly synchronize with the flow of passing messages. The resulting programs enforce the computed order of operations on each resource in the scheduling table, but allow for some timing elasticity: If the execution of an operation takes less than its WCET or WCCT, operations depending on it may start earlier. This elasticity does not compromise the worst-case timing guarantees computed by the mapping tool.

1) Memory handling: Our real-time scheduling and timing analysis use conservative WCET estimations for the (parallel) execution of data-flow blocks on the computing tiles, in isolation. Uncontrolled memory accesses coming from other tiles during execution could introduce supplementary delays that are not taken into account in the WCET figures or by our scheduling tool.

To ensure the timing correctness of our real-time scheduling technique, we need to ensure that memory accesses coming from outside a tile do not interfere with memory accesses due to the execution of code inside the tile. This is done by exploiting the presence of multiple memory banks on each tile. The basic idea is to ensure that incoming DMA transfers never use the same memory banks as the code running at the same time on the CPUs. Of course, once a DMA transfer is completed, the memory banks it has modified can be used by the CPUs, the synchronization being ensured through the use of locks.

We currently ensure this property at code generation time, by explicitly allocating variables to memory banks in such a way as to exclude contentions. While not general, this technique worked well for our case studies, which allowed us to focus the work of this paper on dealing with the NoC resource allocation problem. We are currently working on integrating RAM bank allocation within the mapping algorithm.

VII. EVALUATION

We have evaluated our mapping and code generation method on two applications: The platooning application described in Section V and a parallel Cooley-Tukey implementation of the integer 1D radix 2 FFT over 2 14 samples [START_REF] Bahn | Parallel FFT algorithms on network-on-chips[END_REF]. We chose these two applications because they allow the computation of tight lower bounds on the execution cycle latency and because (for the FFT) an MPPA mapping already exists. This allows for meaningful comparisons, while no tool equivalent to ours was available for evaluation.

For the FFT, we followed the parallelization scheme used in [START_REF] Bahn | Parallel FFT algorithms on network-on-chips[END_REF], with a block size of 2 11 , resulting in a total of 32 tasks. Evaluation is done on the 3x4 MPPA pictured in Fig. 2, where we assume that input data arrives on Tile(0, 0) and the results are output by Tile [START_REF] Amarasinghe | Language and compiler design for streaming applications[END_REF][START_REF]Pnoc: a flexible circuit-switched noc for fpga-based systems[END_REF].

For both applications, after computing the WCET of the tasks and the WCCT of the data transmissions, the mapping tool was applied to build a running implementation and to compute execution cycle latency and throughput guarantees. Then, the code was run, and its performances measured. This allowed us to check the functional correctness of the code and to determine that our tool produces very precise timing guarantees. Indeed, the difference between predicted and observed latency and throughput figures is less than 1% for both examples, which is due to the precision of our mapping algorithms and to the choice of a very predictable execution platform.

The generated off-line schedule (and the resulting code) has good real-time properties. For both the CyCab and the FFT, we have manually computed lower bounds on the execution cycle latency. 4 The lower bounds computed for the CyCab and FFT examples were lower than the latency values computed by our algorithms by respectively 8.9% and 3.4%. For the FFT example, we have also compared the measured latency of our code with that of a classical NoC-based parallel implementation of the FFT [START_REF] Bahn | Parallel FFT algorithms on network-on-chips[END_REF] running on our architecture. For our code, the NoC was statically scheduled, while for the classical implementation it was not. Execution results show that our code had a latency that was 3.82% shorter than the one of the classical parallel FFT code. In other words, our tool produced code that not only has statically-computed hard real-time bounds (which the hand-written code has not) but is also faster.

Our mapping heuristics favor the concentration of all computations and communications in a few tiles, leaving the others free to execute other applications (as opposed to evenly spreading the application tasks over the tiles). The code generated for Cycab has a tile load of 85%-99% for 6 of the 12 tiles of the architecture, while the other tiles are either unused or with very small loads (less than 7%). Using more computing tiles would bring no latency or throughput gains because our application is limited by the input acquisition speed. In the FFT application the synchronization barriers reduce average tile use to 47% on 8 of the 12 MPPA tiles. Note that the remaining free processor and NoC time can be used by other applications.

Finally, we have measured the influence of static scheduling of NoC communications on the application latency, by executing the code generated for Cycab and the FFT with and without NoC programming. For Cycab, not programming the NoC results in a speed loss of 7.41%. For the FFT the figure is 4.62%.

We conclude that our tool produces global static schedules of good quality, which provide timing guarantees close to the optimum.

VIII. CONCLUSION

We have defined a new technique and tool for the offline real-time mapping of data-flow graphs onto many-core architectures supporting programmed arbitration of NoC communications. In doing this, we have shown that taking into account the fine detail of the hardware and software architecture allows off-line mapping of very good precision, even when low-complexity scheduling heuristics are used in order to ensure the scalability of the approach. Our tool synthesizes code with static real-time guarantees that runs faster than (simple) hand-written parallel code.

For the future, our first objective is to perform memory resource allocation at scheduling time (and not during code generation). More generally, our mapping technique could benefit from/to previous work on the scheduling of data-flow specifications and on compilation, but complex evaluation is needed to determine which algorithms scale up to take into account the low-level architectural detail.

Fig. 1 .

 1 Fig. 1. Global flow of the proposed mapping technique

Fig. 2 .

 2 Fig. 2. A 4x3 tile MPPA platform with 2D mesh NoC. Black rectanges are the NoC routers.

Fig. 3 .

 3 Fig. 3. Generic router for a 2D mesh NoC with X-first routing policy

Fig. 4 .

 4 Fig. 4. The computing tile of our MPPA architecture

Fig. 5 .

 5 Fig. 5. Data-flow graph of a platooning application

Fig. 6 .

 6 Fig. 6. Scheduling table covering one communication path on our NoC. Only the 6 resources of interest are represented (out of 70). Time flows from top to bottom.

  DataSize : worst-case size of data to transmit FreeIntervalList : list of free intervals before reservation ScheduleLength : schedule length before reservation Output: ReservedIntervalList : reserved intervals NewIntervalList : list of free intervals after reservation NewScheduleLength : schedule length after reservation 1: NewIntervalList := FreeIntervalList 2: ReservedIntervalList := ∅ 3: while DataSize > 0 and NewIntervalList = ∅ do 4:

  Path : list of resources (the communication path) StartDate : date after which the data can be sent DataSize : worst-case data size (in 32-bit words) Input/Output: SchedulingTable : scheduling table 1: for i := 1 to length(Path) do

	Procedure 1 MapCommunicationOnPath
	2:	ShiftSize := (i -1) * SegmentBufferSize;
	3:	FreeIntervalList[i] :=
		GetIntervalList(SchedulingTable, GetSegment(Path,i), ShiftSize)
	4:	ShiftedIntervalList[i] :=
		ShiftLeftIntervals(FreeIntervalList[i],ShiftSize)
	5: end for
	6: PathFreeIntervalList :=IntersectIntervals(ShiftedIntervalList);
	7: (ReservedIntervals,NewIntervalList,NewScheduleLength) :=
		ReserveIntervals(DataSize,PathFreeIntervalList,
		length(SchedulingTable));
	8: (IntervalForLock,NewIntervalList,NewScheduleLength) :=
		ReserveIntervals(LockPacketLength,NewIntervalList,
		NewScheduleLength);
	9: ReservedIntervals := AppendList(ReservedIntervals,IntervalForLock)
	10: for i := 1 to length(Path) do
	11:	ShiftSize := (i-1)*SegmentBufferSize;
	12:	FinalIntervals[i] :=ShiftRightIntervals(ReservedIntervals,ShiftSize);
	13: end for
	14: if NewScheduleLength > length(SchedulingTable) then
	15:	SchedulingTable :=
		IncreaseLength(SchedulingTable,NewScheduleLength);
	16: end if
	17: SchedulingTable :=
		UpdateSchedulingTable(SchedulingTable,Path,FinalIntervals

Input:

Implementing an optimal control scheme for such an application may require more resources than the application itself, which is why dynamic/online scheduling techniques are often preferred.

Throughput in this context means the number of execution cycles started per time unit. It can be different from the inverse of the latency because we allow one cycle to start before the end of previous ones, provided that data-flow dependencies are satisfied.

To compute these lower bounds we simplify the hardware model by assuming that the resources N (i, j)(k, l) generate no contention (i.e. they allow the simultaneous transmission of all packets that demand it). We only take into account the sequencing of operations on processors and DMAs and the contentions on resources In(i, j).


Frédéric Giroire

Stéphane Pérennes

Issam Tahiri

Grid spanners with low forwarding index for energy efficient networks *

Keywords: spanning subgraphs, forwarding index, energy saving, routing, grid sous-graphes couvrants, indice de transmission, efficacité énergétique, routage, grille

. Motivated by energy efficiency, we look, for different numbers of edges, at the best spanning graphs of a square grid, namely those with a low forwarding index.

Introduction

A routing R of a given connected graph G of order N is a collection of N (N -1) simple paths connecting every ordered pair of vertices of G. The routing R induces on every edge e a load that is the number of paths going through e. The edge-forwarding index (or simply the forwarding index) π(G, R) of G with respect to R is the maximum number of paths in R passing through any edge of G. It corresponds to the maximum load over all edges of the graph when R is used. Therefore, it is important to find routings minimizing this index. The forwarding index π(G) of G is the minimum π(G, R) over all routings R's of G. This parameter has been studied for different graph classes (examples can be found in [START_REF] Bouabdallah | On the edge forwarding index problem for small graphs[END_REF], [START_REF] Fernandez De La Vega | The forwarding indices of random graphs[END_REF], [START_REF] De La Vega | The forwarding index of communication networks with given connectivity[END_REF]) and this survey [START_REF] Xu | The forwarding indices of graphs -a survey[END_REF] gives a global view on the known results.

We call a connected spanning subgraph of a graph G, a spanner of G. More precisely, it is a connected subgraph that has the same set of vertices as G. Our goal is to find, for a given bound on the number of edges, the best spanner of G, namely the one with the minimum forwarding index. The problem can also be viewed as: for a given bound U on the forwarding index, find a spanner F of G with minimum number of edges such that π(F ) ≤ U .

Knowing how to solve this problem is very interesting in practice for network operators willing to reduce the energy consumed by their networks. In fact, most of the network links consume a constant energy independently of the amount of traffic they are flowing [START_REF] Chabarek | Power awareness in network design and routing[END_REF], [START_REF] Restrepo | Energy profile aware routing[END_REF]. Therefore, it was proposed to reduce the energy used by the network links by turning some of them off, or more conveniently, putting them into an idle mode. Outside the rush hours, several studies [START_REF] Chiaraviglio | Energy-aware umts core network design[END_REF], [START_REF] Chiaraviglio | Reducing power consumption in backbone networks[END_REF], [START_REF] Baliga | Energy consumption in IP networks[END_REF], [START_REF] Giroire | Minimizing routing energy consumption: from theoretical to practical results[END_REF], [START_REF] Giroire | Minimization of network power consumption with redundancy elimination[END_REF], [START_REF] Giroire | Optimizing rule placement in software-defined networks for energy-aware routing[END_REF], show that a good choice of the links to turn off can lead to significant energy savings, while keeping the same communication quality. In the case where the throughputs from every node to every other node are of the same order, and where the capacities also lie in the same small range, a good choice of those links is reduced to the problem of finding spanners of the network with low forwarding indices.

In this paper, we consider the case in which the initial graph is a square grid. We consider the asymptotic case with n large. We have two main contributions.

On one side, we know that the forwarding index of the n × n grid G n is n 3 2 , see Proposition 1 [START_REF] Giroire | Energy Efficient Routing by Switching-Off Network Interfaces[END_REF]. G n has 2(n -1) 2 ∼ 2n 2 edges. An important remark is that the load of the edges is lower in the corner than in the middle of the grid. Using that, we show that we can build spanners of G n with much fewer edges (only 13/18 ≈ 72% of the edges) and the same forwarding indices as G n . We show that the proposed spanners are close to optimum in the sense that we prove that it is impossible to build spanners with fewer than 4/3n 2 edges (66% of the edges).

On the other side, the smallest possible spanner of the n × n grid G n is a spanning tree. The forwarding index of the best spanning tree is asymptotically 3n 4 8 , see Proposition 2 [START_REF] Giroire | Energy Efficient Routing by Switching-Off Network Interfaces[END_REF]. When we add edges and consider spanners with a larger number of edges, the load on the edges decreases, and so does the forwarding index. In this paper, we study how the forwarding index decreases, when we increase the number of edges.

The following table summarizes our results:

Spanning tree Spanners Grid forwarding index Proposition 1 [START_REF] Giroire | Energy Efficient Routing by Switching-Off Network Interfaces[END_REF] The forwarding index of G n is asymptotically n 3 2 .

3 8 n 4 1 2a n 4 (2 ≤ a < n) 1 2 n 3 1 2 n 3 lower
Proposition 2 [START_REF] Giroire | Energy Efficient Routing by Switching-Off Network Interfaces[END_REF] For n ≥ 3, The spanning tree of G n with the minimum forwarding index is a tree with centroid of degree 4 and 4 branches of almost equal sizes. its forwarding index is asymptotically 3n 4 8 .

2 Spanners with the forwarding index of the grid, n 3 2 , but much fewer edges

In this section, we first show that a spanner with the forwarding index of the grid has at least 4n 2 3 = 12n 2 9 edges. We then provide spanners with 13n 2 9 edges. But, before, we present some notations that will be used throughout the paper.

Notations We note by G n = (V n , E n ) the n × n square grid, where V n is the set of Vertices and E n is the set of edges. A square grid can always be seen as n rows intersecting n columns. We name v(r, c) the vertex at the intersection of row r ∈ [n] with column c ∈ [n], where [n] denotes the interval of the integer numbers between 1 and n. An edge joining v(r, c) to v(r, c + 1) is named e h (r, c) and an edge joining

v(r, c) to v(r + 1, c) is named e v (r, c). Proposition 3 For any F spanner of G n such that π(F ) ≤ n 3
2 , F must have, asymptotically, at least 4n 2 3 edges.

Proof: Consider F a spanner of G n and let R be a routing of F such that π(F, R) ≤ n 3 2 . For an integer l ∈ [n], we call load on line l, the sum of the load on the edges e v (l, j) ∈ E(F ), for j ∈ [n]. The load on line l is 2l(n -l)n 2 as there are ln vertices over line l and (n -l)n vertices below. If F has n -x l edges on line l, there exists at least one of these edges with load at least 2l

(n-l)n 2 n-x l . As π(F, R) ≤ n 3 2 , we should have 2l(n -l)n 2 n -x l ≤ n 3 2 .
That is

n -x l ≥ 4l(n -l) n .
Thus, F should have at least

n l=1 4l(n-l) n
vertical edges. The same argument independently holds for the horizontal edges. Hence, a spanner of the grid, with load lower than n 3 2 on all edges, has at least

2 n l=1 4l(n -l) n ≈ 4n 2 3 edges.
Theorem 1 There exists F n a spanner of G n such that π(F n ) ∼ n 3 2 and its number of edges is asympotically equal to 13n 2 9 . Proof: Let us first explain the intuition behind the construction of the spanner of the grid, F n . We know from the proof of Proposition 3 the ratio of edges needed in every row or column in order to satisfy the lower bound. We cut the grid into small squares. Then, according to the position of the square, we use only the number of needed horizontal edges and vertical edges in each square according to the lower bound. It turns out that adding only few edges to ensure the connectivity is enough to get a spanner F n with a routing R such that π(F n , R) ∼ n 3 2 . Construction of F n . Let k be an integer number such that 1 ≤ k ≤ n. We divide G n into small square grids of size k × k. We do so by partitioning vertices of

G n into ( n k ) 2 sets S (i,j) with i ∈ [ n k ] and j ∈ [ n k ]: S (i,j) = {v(r, c) ∈ V n ; i -1 < r k ≤ i, j -1 < c k ≤ j}.
We call a vertex in S (i,j) that has a neighbour in G n outside S (i,j) a border vertex. Let us now describe a spanner F n that verifies our theorem. An example of it is shown in Figure 1 in the case of n = k 2 = 7 2 . Let t be the function defined on integers by t(x) = 4xk(n -xk)k/n 2 . It represents the number of needed columns (respectively rows) for a square that is on the x-th position horizontally (respectively vertically). We build F n starting from a subgraph that has all vertices of G n and no edges. For every S (i,j) , i, j ∈ [ n k ], we choose edges to connect vertices in S (i,j) in the following way:

-we add to F n all edges e v (r, c) such that (r mod k) ∈ {1, . . . , t(i)} (red edges in Figure 1) and -all edges e h (r, c) such that (c mod k) ∈ {1, . . . , t(j)} (blue edges in Figure 1); -then we add to F n simple paths just to connect the remaining independent vertices (green edges in Figure 1).

-We then add all edges that do not have both endpoints in the same set S (i,j) (black edges in 1). We show in the following that adding all of them is not strictly necessary.

Description of the routing R. We now give a routing of the spanner F n , R. For every ordered pair of vertices (v(r a , c a ), v(r b , c b )) of V n , we describe the path connecting v(r a , c a ) to v(r b , c b ) in R. We distinguish two types of ordered pairs of vertices:

• Type-1 pairs: r a /k = r b /k or c a /k = c b /k . Notice that this type includes ordered pairs with vertices that belongs to the same set S (i,j) .

• Type-2 pairs: All the ordered pairs that do not belong to the first type.

For the Type-1 pairs, R uses the shortest path routing. For Type-2 pairs, R uses a three-segment path.

An example of such path is shown in Figure 1. We name

i a = r a /k , i b = r b /k , j a = c a /k and j b = c b /k : • Step-1: Let i m = min(i a , i b , n/k -i a , n/k -i b ) and j m = min(j a , j b , n/k -j a , n/k -j b ).
The first segment is the shortest path from v(r a , c a ) to one of the two border vertices of S (ia,ja) that are on row k(i a -1) + t(j m ). Among the two vertices, we choose v(r x , c x ), which has the smallest distance to S (ia,j b ) (as the first black vertex on the route in Figure 1).

• Step-2: Similarly, two border vertices of S (i b ,j b ) are on column k(i b -1) + t(i m ). Among these two vertices, v(r y , c y ) is the one that has the smallest distance to S (ia,j b ) (as the third black vertex on the route in Figure 1). The second segment will be linking v(r x , c x ) to v(r y , c y ) by using the path [v(r x , c x )v(r x , c y )v(r y -c y )], which is the shortest path from v(r x , c x ) to v(r x , c y ) composed of the two direct paths [v(r x , c x )v(r x , c y )], following row r x , and [v(r x , c y )v(r y -c y )], following column c y .

• Step-3: The third and last segment will be the shortest path from v(r y , c y ) to v(r b , c b ).

Note that k may be an arbitrary integer between 1 and n. We choose a k such that 1 k √ n. For instance, we may choose k = n 1/3 . Number of edges of F n . Let us compute the number of edges in the spanner, F n . First the edges used in the subgraph induced by S (i,j) are all the edges on a row from 1 to t(i), all edges on a column from 1 to t(j), to which we add the edges that connect the rest of vertices through a spanning tree. Hence the number of edges in S (i,j) is:

≈ t(i) • k + t(j) • k + (k -t(i))(k -t(j)) ≈ k 2 + t(i)t(j) (1)
≈ k 2 (1 + 16ijk 2 (n -ik)(n -jk) n 4 ) ≈ k 2 (1 + 16ijk 2 n 2 + 16i 2 j 2 k 4 n 4 - 16i 2 jk 3 n 3 - 16ij 2 k 3 n 3 )

Inria

The sum of those edges considering all the subsets S (i,j)

(with i ∈ [ n k ] and j ∈ [ n k ]) is : ≈ k 2 n/k i=1 n/k j=1 (1 + 16ijk 2 n 2 + 16i 2 j 2 k 4 n 4 - 16i 2 jk 3 n 3 - 16ij 2 k 3 n 3 ) ≈ k 2 [ n 2 k 2 + 16k 2 n 2 ( i i) 2 + 16k 4 n 4 ( i i 2 ) 2 -2 • 16k 3 n 3 ( i i 2 )( i i)] ≈ k 2 [ n 2 k 2 + 16k 2 n 2 • n 4 4k 4 + 16k 4 n 4 • n 6 9k 6 -2 • 16k 3 n 3 • n 5 6k 5 ] ≈ k 2 [ n 2 k 2 + 4n 2 k 2 + 16n 2 9k 2 - 32n 2 6k 2 ] = 13 9 n 2 + o(n 2 )
The number of the remaining edges is ≈ 2

n 2 k = o(n 2 ), as k 1.
Therefore, as stated in the theorem, the number edges of F n is asymptotically equal to 13n 2 9 + o(n 2 ).

Load of the edges of F n . Lets now verify that every edge has an asymptotic load which is not greater than n 3 2 + o(n 3 ). Consider an edge e h (r, c) whose incident points are in S (i,j) . The number of Type-1 pairs that may use e h (r, c) is bounded by the number of pairs having one endpoint in S (i,ja) and S (i,j b ) for some j a , j b ∈ [ n k ] and those having one end point in S (ia,j) and

S (i b ,j) for some i a , i b ∈ [ n k ]. The number of these pairs is bounded by 2k 2 n 2 = o(n 3 ) (as k 2 = o(n)).
Then, for Type-2 pairs, we can start by the load induced by the segments of paths described previously in step-1 and step-3. This load is clearly bounded by the number of pairs having one endpoint inside S (i,j) and another endpoint outside S (i,j) . The number of these pairs is bounded by: 2k 2 (n -k) 2 = o(n 3 ) (as k 2 = o(n)). For Step-2, as the construction of the spanner F n has the needed density of edges, the load is kept below n 3 2 + o(n 3 ). The same argument holds for the black edges between two adjacent subsets S (ia,ja) and S (i b ,j b ) . This ends the proof.

3 Spanners with forwarding indices in the range ] n 3 2 , 3n 4

[ and Lower bounds

We first provide spanners with forwarding indices in the range ] n 3 2 , 3n 4 8 [ in Proposition 4. We then prove that these spanners have a number of edges of the optimum order, see Proposition 5

Spanners' constructions

Proposition 4 Let a be an integer such that, 2 ≤ a ≤ n. There exists a spanner F n (a) of G n with asymptotically n 2 + 4 9 a 2 edges and π(F n (a)) ≤ n 4 2a .

Proof: We build a spanner of G n , F n (a), in the following way. We divide the grid into a 2 sectors. A point is in Sector (i, j) if its coordonnates in the grid (x,y) are such that n a i ≤ x < n a (i + 1) and n a j ≤ y < n a (j + 1). Each of these sectors has (n/a) 2 vertices. We call center of the sector (i, j) the vertex ((i + 1/2)n/a, (j + 1/2)n/a). We consider the a × a subgrid linking all the sectors' centers. We then connect all the remaining vertices of a sector to its center with a spanning tree. This way, we get F n (a). We now build a routing R for F n (a). The demand between two vertices of the same sector are routed on the tree spanning their sector using the unique shortest path between them. The demand between two vertices of different sector is first routed to their centers, and then is routed in the a × a grid.

Let us compute the load of the routing R. We first consider the edges of the a × a subgrid. We know that an a × a grid has a routing with load a 3 /2 (Proposition 1). Thus, we know that it also has a w-routing of load wa 3 /2. Each vertex of the a × a grid receives the load of the (n/a) 2 vertices connected to it. Thus, we take w = (n/a) 2 and we obtain a w-routing of the a × a grid of load a 3 2 ( n a ) 4 = n 4 2a . We then consider an edge that does not belong to the a × a grid. The only paths that can use this edge are paths going from any vertex of the grid to a vertex of its sector. Thus, its load is smaller that (n/a) 2 n 2 = n 3 a 2 . This load is smaller than the maximum load on the a × a grid as soon as a 2 ≥ 2a which means as soon as a ≥ 2.

Therefore π(F n (a), R) = n 4 2a . Let us now consider the number of edges of the spanner F n (a). The number of edges necessary to connect all the nodes is n 2 -1. If we choose well these edges, we just have to add a 2 edges to obtain the a × a grid (see Figure 2, additional edges are in red). F n (a) thus has n 2 + a 2 edges. We can improve the Inria

Conclusion

We succeeded at providing spanners of the n×n grid with a small number of edges for a given forwarding index. Such spanners are important for energy efficient networks in which the traffic has to be routed in the network while using a minimum number of equipments. The unused equipments are then turned off to save energy. We leave as open two problems.

We propose spanners with a number of edges of optimum order for a forwarding index. More precisely, we have provided spanners of the grid with n 2 + 4 9 a 2 edges and forwarding indices 1 2a n 4 (2 ≤ a < n). We proved that it is impossible to have spanners with the same FI and fewer than n 2 + 4 9 (0.1a) 2 edges. It would be very nice to succeed in filling the gap between the lower bounds and the constructions.

Similarly, we describe spanners with 13/9n 2 edges and with the same forwarding index of the full grid G n . We proved that spanners with such a forwarding index should have at least 12/9n 2 edges. Would it be possible to find spanners with such a number of edges?

Last, we focused on a specific network in this work, the square grid. We are also interested by more general graphs. In particular, the arguments to derive lower bounds can be used for more general planar graphs with bounded degrees. It would be very interesting to find results and constructions for other families of planar graphs.

Figure 1 :

 1 Figure 1: Construction of the spanner F n of Theorem 1, for n = 7 2 , and an example of path of the routing R of F n (from the yellow vertex to the pink vertex).

  Figure 2 provide a sketch of the construction of the spanner.

Figure 2 :

 2 Figure 2: Spanner of Proposition 4 for n = 21 and a = 3. Edges of the a × a grid are in bold. Edges that are not in a spanning tree of G n are in red. Sectors with (n/a) 2 = 7 2 vertices are separated by dashed gray lines.

RR n°8643

* This work has been partially supported by ANR project Stint under reference ANR-13-BS02-0007, ANR program Investments for the Future under reference ANR-11-LABX-0031-01, ANR VISE, ANR Skyflow, CNRS-FUNCAP project GAIATO, the associated Inria team AlDyNet, the project ECOS-Sud Chile and Paca Region.

spanner by using the results of Section 2. In Theorem 1, we show that we can find a spanner of an a × a subgrid with 13 9 a 2 edges and a routing R with the same load as a full grid with 2a 2 edges. By doing so, we get a new spanner F n (a), with n 2 + 4 9 a 2 edges and π(F n , R ) = n 4 2a .

We can rewrite the result of Proposition 4 to point out the impact of additional edges in general (Corollary 1) and when we start from a spanning tree (Corollary 2).

Corollary 1 There exist:

-A spanner of G n with n 2 + p 2 edges, and an asymptotic forwarding index of 8 , that is a factor α less than the one of the optimum spanning tree, while using 64 81 α 2 0.79α 2 additional edges compared to a spanning tree. Proof: Recall that an optimum spanning tree has forwarding index 3n 4 8 , see Proposition 2. Dividing it by α means getting the forwarding index 3n 4 8α = n 4 2(4α/3) . This is achieved by the spanner F n (a), with a = 4α/3. The spanner has an additional number of edges compared to the spanning tree equal to 4 9 (4α/3) 2 0.79α 2 .

Lower bounds

Proposition 5 There exist no spanners of G n with n 2 + p 2 edges and a forwarding index less than Proof: Let us consider a spanner of G n that has n 2 + p 2 edges. We build a multigraph in the following way. We start by assigning to every node a weight of 1. Then, while there is still a vertex with degree 1 or 2, we delete this vertex and the edges connecting it to the graph and divide its weight evenly among its neighbors; in case the removed vertex was of degree 2, we also connect the two neighbors afterwards. At the end of this process, we get a multigraph H such that the number of its vertices N and the number of its edges M are related by the following equation: N + p 2 = M . Indeed, every time a vertex is removed in the process leading to H, the number of edges is decreased by 1. Since all the vertices in H have a degree strictly greater than 2, we have 3 2 N ≤ M . This implies with the previous equation that

Hence, we have N ≤ 2p 2 . Notice that the total weight is equal to n 2 . We now apply the weighted version of the planar separator theorem [START_REF] Djidjev | On the problem of partitioning planar graphs[END_REF] on H: there exists a partition of the vertices of H into three subsets A, S, and B, such that each of A and B has at most a weight 2n 2 /3, S has less than √ 6 2p 2 vertices (The original graph is of a bounded degree 4.) and there are no edges with one endpoint in A and another endpoint in B. This directly gives an edge cut of the original graph which has less than 2 √ 6 2p 2 edges and which partitions the original graph's vertices into two subsets of size at most 2n 2 /3. Therefore, any routing of this spanner will induce, at least on one edge of the cut, a load that is greater than:


OptiDis: A parallel Fast Multipole Dislocation Dynamics code

Pierre Blanchard 1 -Arnaud Etcheverry 1 -Olivier Coulaud 1 -Laurent Dupuy 2 -Marc Blétry 2 -Eric Darve 3

[1]: HiePACS team project, Inria Bordeaux -Sud-Ouest, 200, rue Vieille Tour 33405 Talence Cedex France | [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF]: Service de Recherches Métallurgiques Appliquées, CEA-Saclay, 91191 Gif-sur-Yvette, France | [START_REF] Fong | The black-box fast multipole method[END_REF]: Mechanical Engineering, Stanford University, CA, USA ABSTRACT Among all the steps involved in DD simulations, the computation of the internal elastic forces and energy are by far the most resources consuming. However, since these are long-ranged and fast decreasing interactions, hierarchical algorithms like the Fast Multipole Method (FMM) are well suited for their fast evaluation. The relatively low accuracy required for the interaction forces between dislocations brought us to develop a more efficient approximation method for the farfield. On the other hand, the nearfield interactions are still evaluated analytically, which required a rather performant implementation (AVX, GPU...).

Regarding parallelism, our code benefits from a hybrid OpenMP/MPI paradigm and a cache-conscious datastructure. Finally, an accurate handling of topological elements intersecting the structure of the octree was considered. The latter feature implied careful modifications of the P2M/L2P operators in order to deal with shared memory model of parallelism.

MODEL

The motion of dislocations is ruled by a prescribed mobility law (e.g. viscous glide) depending on the internal (i.e. created by dislocations themselves) and external (i.e. applied by user) contributions on the nodal forces. The isotropic elastic stress field created by a dislocation loop (b , t ) at point x in space is given by Mura's formula [START_REF] Mura | Micro-mechanics of Defects in Solids[END_REF]

σ ij (x) = µ 8π ((A ij )(x) + (A ji )(x) + 2 1 -ν (B ij -δ ij B pp )(x))
where (µ, ν) are the Lamé coefficients and

(A ij )(x) = (C ) R ,ppm (x, y)ε jmk b k t i dy (B ij )(x) = (C ) R ,ijm (x, y)ε nmk b k t n dy with R(x, y) = |x -y|.
The nodal force f e n acting at the extremities of a finite dislocation line (b, t) is obtained by integration of the Peach-Koehler force (f

P K = (σ(x) • b) × t(x)) over the target line, i.e. (f e n ) α = (C) ε αβγ σ βp b p t γ N n (x)dx
where N n=1,2 are linear shape functions. While the isotropic elastic interaction energy is given by

E(C , C) = - µ 8π (C) (C ) R ,kk (b i dx i b j dx j - 2µ 1 -ν (b i dx j b i dx j -b j dx i b i dx j ))+R ,ij 2µ 1 -ν b j dx k b i dx k
The cost of updating the nodal forces (and the energy) at each time step is quadratic and involves the evaluation of 2 line integrals, therefore it is usually the bottleneck of DD simulations.

A NEW FFT-ACCELERATED FMM

Fast Multipole DD

The idea of the FMM is to balance the computational cost of near and farfield interactions by approximating the farfield and rely on a tree structure (see figure 1) to perform all computations. In our case we consider the interaction potential R(x, y) = {R ,ij (x, y)} i,j=1..3 and apply the interpolation formula

R(x, y) ≈ |α|≤p S α (x) |β|≤p R(x α , ȳβ )S β (y)
where xα and ȳβ denote 2 sets of interpolation nodes in 2 well separated clusters containing x and y. S denotes a polynomial interpolator (either Lagrange or Chebyshev). The farfield approximation consists in 3 steps: P 2P Analytical evaluation of the nearfield.

P 2P ↑ M 2M/L2L ↓ M 2L ↑ P 2M/L2P ↓
P 2M/L2P Integration of polynomial interpolators on all leaf cells returns expansion:

M β = (C)
S β (y)dy M 2M/L2L Transfer expansions between levels. M 2L Transfer expansions between cells in interaction by applying

L α = R αβ S β
Both stress and energy expansions share the same M2L operators R αβ , but the extra derivation involved in the stress field is applied to the interpolator. This results in a slightly lower accuracy but less M2L operations to be performed (see figure 2). interpolation order: p relative error

|F | using R ,ijm |F | using R ,ij ∇ m W using R ,ij

Acceleration by Fast Fourier Transform (FFT)

Even though the Chebyshev-based FMM [START_REF] Fong | The black-box fast multipole method[END_REF] provides a very accurate algorithm it is very expensive in terms of memory and CPU time, namely O(p 6 ) for the M2L step. Thus a new interpolation scheme based on an equispaced grid (Lagrange) was considered. In this approach R αβ is Toeplitz since R = R(|x -y|).

As illustrated below, R αβ can be seen as a convolution and applied in Fourier space with a linear cost. The scheme remains stable and sufficiently accurate in the scope of DD simulations (min. accuracy on force computation ≈ 10 -3 ), while dramatically decreasing the computational cost of the M2L, namely to O((2p -1) 3 ). Let p = 5 then in double precision the memory required for the M2L equals:

316 × (2 × 5 -1) 3 × 7 × (2 × 8) ≈ 25 MBytes

IMPLEMENTATION AND PERFORMANCES

Our experimentations were performed on the core program OptiDis whose data structure relies heavily on the open source ScalFMM library [START_REF]Scalfmm: software library to simulate large scale n-body interactions using the fast multipole method[END_REF]. The latter also provides the generic Fast Multipole algorithms. Op-tiDis implements its own dislocation specific P2M/L2P routines involving exact integration over segments as well as analytic expressions for nearfield taken from Arsenlis et al [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF].

Figure 1 : Plop. 1 Figure 1 :

 111 Figure 1: Plop.

Figure 1 :Figure 2 :

 12 Figure 1: Accuracy of the nodal force evaluation for different interpolation schemes with respect to the interpolation order. The topological configurations are illustrated figure ??.

  RX=FFT -1 (FFT(C 0 ) : FFT(X)) with C 0 = a b c e d

1 Figure 3 :

 13 Figure 3: Near-and Farfield computational time balancing for a uniform distribution of dislocation loops and increasing tree depth (left). An example of defects distribution and a propagating Frank-Read source (right).

Figure 4 :

 4 Figure 4: Evolution of 6 screw Frank-Read sources (top) inside a grain of Zr with a large density of defects (d = 1.10 22 m -3 ). Elementary events (bottom): (a) double jog (b) dissociation.


Pierre Blanchard

Arnaud Etcheverry

Olivier Coulaud

Laurent Dupuy

Marc Bletry

Marc Bletry Optidis

OptiDis: Toward fast anisotropic dislocation dynamics based on Stroh formalism

come   L'archive ouverte pluridisciplinaire

OptiDis: Toward fast anisotropic DD based on Stroh formalism.

Pierre Blanchard 1 -Arnaud Etcheverry 1 -Olivier Coulaud 1 -Laurent Dupuy 2 -Marc Blétry 2

[1]: HiePACS team project, Inria Bordeaux -Sud-Ouest, 200, rue Vieille Tour 33405 Talence Cedex France | [START_REF] Aubry | Use of spherical harmonics for dislocation dynamics in anisotropic elastic media[END_REF]: Service de Recherches Métallurgiques Appliquées, CEA-Saclay, 91191 Gif-sur-Yvette, France ABSTRACT Dislocation Dynamics (DD) simulations in the hypothesis of isotropic elasticity have proved great reliability in order to predict the plastic behaviour of crystalline materials. However it is often the case at high temperature (for instance in irradiated BCC iron) that the structural properties of a material will be better described using full anisotropic treatment of the elastic interaction between dislocations. The computation of the internal elastic forces is by far the most resources consuming step in DD simulations, which is even more true for anisotropic elasticity in the absence of explicit Green's function.

L. Dupuy, J. Soulacroix and M. Fivel showed that the approaches summarized in Yin [START_REF] Yin | Efficient computation of forces on dislocation segments in anisotropic elasticity[END_REF] can be accelerated using spherical harmonics expansions of the Stroh matrices. This feature was implemented in the DD code OptiDis in order to power the anisotropic forces computation. Here we recall the formalism and we discuss optimizations, performances as well as motivations for future developments.

ANISOTROPIC MODEL

The stress field created by a dislocation loop (b ′ , t ′ ) at field point x is given by Mura's formula [START_REF] Mura | Micro-mechanics of Defects in Solids[END_REF]

σ js (x) = ε ngr C jsvg C pqwn b ′ w C ∂G vp ∂x q (x -x ′ )dx ′ r (1)
The nodal force f e n acting at the extremities of a finite dislocation line (b, t) is obtained by integration of the Peach-Koehler force (f

P K = (σ(x) • b) × t(x)) over the target line, i.e. (f e n ) α = (C) ε αβγ σ βp b p t γ N n (x)dx
where N n=1,2 are linear shape functions. The cost of updating the nodal forces at each time step is quadratic and involves the evaluation of 2 line integrals, therefore it is usually the bottleneck of DD simulations. Moreover there does not exist an analytic closed form for the anisotropic elastic Green's function G. Recently Aubry et al. [START_REF] Aubry | Use of spherical harmonics for dislocation dynamics in anisotropic elastic media[END_REF] used an integral representation of G and developped a fast method based on spherical harmonic expansions in order to evaluate the double line integral semi-analytically. On the other hand, past works [START_REF] Yin | Efficient computation of forces on dislocation segments in anisotropic elasticity[END_REF] showed that the anisotropic stress field can be efficiently described using the Stroh sextic formalism combined with (2) the Willis-Steeds-Lothe formula for the finite line. In the collinear case (d = 0) the expression is slightly more complicated but can be condensed as follows

Willis-Steeds-Lothe

u m,s = 1 4πd ε jsn b i C ijkl t n -m l Q mk + n l (nn) -1 (nm)Q mk + n l (nn) -1 S T mk τ 2 τ 1 (2)
u m,s (d = 0) = 1 r 2 - 1 r 1 1 4π b i Z ims
The singularity in the limit r → 0 is currently handled using a simple cutoff parameter like the one defined in [START_REF] Yin | Efficient computation of forces on dislocation segments in anisotropic elasticity[END_REF].

Anisotropy ratio

The degree of anisotropy is quantified by the ratio A = 2C 44 /(C 11 -C 12 ). For the BCC α -F e, this ratio goes from A 0K = 2.3 to A 1200K = 7.1.

SPHERICAL HARMONIC ANALYSIS

Q S Z

Spherical harmonic analysis

Stroh matrices only depend on the orientation of the source, i.e. X = X(θ, φ) (see fig 2) hence they can be expanded into spherical harmonics.

X(θ, φ) ≈ ℓ max ℓ=0 ℓ m=-ℓ x m ℓ Y m ℓ (θ, φ)
where Y m ℓ denotes the well known spherical harmonics and Given that Stroh matrices are real valued the expansion reduces to

x m ℓ = π 0 2π 0 X(θ, φ)Y m ℓ (θ, φ) sin θdθdφ (
X ≈ ℓ max ℓ=0 x 0 ℓ Y 0 ℓ + 2 ℓ m=1 Re(x m ℓ )Re(Y m ℓ ) -Im(x m ℓ )Im(Y m ℓ )
On the other hand depending on the symmetries of X in θ or φ some coefficients of the expansions are known to be null (potentially a lot). Once implemented these simplifications lead to a significant acceleration of the method (see fig 4). The precomputation of the 2 integrals (3) is done exactly following the method described by Driscoll and Healy [START_REF] James | Computing fourier transforms and convolutions on the 2-sphere[END_REF], namely using an equispaced quadrature in φ and a gaussian quadrature in cos θ.

IMPLEMENTATION AND PERFORMANCES

Our experimentations were performed on the core program OptiDis whose data structure relies heavily on the open source ScalFMM library [START_REF]Scalfmm: software library to simulate large scale n-body interactions using the fast multipole method[END_REF]. The latter also provides the generic Fast Multipole algorithms.

OptiDis is a parallel version of NumoDis, it implements almost all functionalities of NumoDis while providing a hybrid OpenMP/MPI paradigm and a cache-conscious data structure. ONGOING & PERSPECTIVES

Ongoing

• Optimized expansion for hexagonal crystallographies

Perspectives

• Implementation of the farfield (either iso-or anisotropic)

• Efficient analytic integration of the expansion over the target segments

• Derivation of a consistent non-singular theory for the Stroh approach

Figure 1 :

 1 Figure 1: Notations for the stress created by segment AB on a field point P .

Figure 2 :

 2 Figure 2: Stroh matrices components at (T = 25 • C) for θ ∈ [0, π] and φ ∈ [0, 2π].

3 )Figure 3 :

 33 Figure 3: Accuracy of the expansion for various temperatures.

Figure 4 :

 4 Figure 4: Theoretical number of flops involved in self force evaluation (left) and relative computational cost for ≈ 3000 segments based on CPU time (right).

Figure 5 :

 5 Figure 5: Isotropic near-and farfield computational time balancing for a uniform distribution of dislocation loops and increasing tree depth (left). An example of defects distribution and a propagating Frank-Read source (right).