Emilien Kofman

A formal approach to the mapping of tasks on an heterogenous multicore, energy-aware architecture

The search for optimal mapping of application (tasks) onto processor architecture (resources) is always an acute issue, as new types of heterogeneous multicore architectures are being proposed constantly. The physical allocation and temporal scheduling can be attempted at a number of levels, from abstract mathematical models and operational research solvers, to practical simulation and run-time emulation.

This work belongs to the first category. As often in the embedded domain we take as optimality metrics a combination of power consumption (to be minimized) and performance (to be maintained). One specificity is that we consider a dedicated architecture, namely the big.LITTLE ARM-based platform style that is found in recent Android smartphones. So now tasks can be executed either on fast, energy-costly cores, or slower energy-sober ones. The problem is even more complex since each processor may switch its running frequency, which is a natural trade-off between performance and power consumption. We consider also energy bonus when a full block (big or LITTLE) can be powered down. This dictates in the end a specific set of requirements and constraints, expressed with equations and inequations of a certain size, which must be fed to an appropriate solver (SMT solver in our case). Our original aim was (and still is) to consider whether these techniques would scale up in this case. We conducted experiments on several examples, and we describe more thoroughly a task graph application based on the tiled Cholesky decomposition algorithm, for its relevant size complexity. We comment on our findings and the modeling issues involved.

I. INTRODUCTION

The multiprocessor scheduling problem is a vast and still open topic. Diverse refinements of this problem exist, and most of them are NP-Hard. However the many and recent improvements in SAT and SMT solvers allow to solve medium size instances of such complexity especially since realistic problems rarely reach the worst case complexity. The demand for more computing power in consumer electronics raises again with the new high definition applications. Increasing the frequency in order to run sequential implementations such as video decoding at an acceptable frame rate is no more a valid option. Using parallel implementations and offloading tasks to low-power resources allows to keep our devices warm and smooth. In this work we show a practical approach of multiprocessor scheduling with an SMT solver under realistic hypothesis. This allows us to actually compile and run an application that benefits from this scheduling and check with This work was partly funded by the French Government (National Research Agency, ANR) through the "Investments for the Future" Program reference #ANR-11-LABX-0031-01. this implementation whether our modeling assumptions are correct or not.

Not every program is a suitable candidate for this methodology. Section III motivates the approach and explains the hypothesis that we need to validate. We then explain the energy model that is used in this approach: given a simple CPU energy model, the instant power of a CPU raises linearly with frequency and with the square of the voltage. Even if those non-linearities introduce yet more practical complexity into the already NP-Hard scheduling problem, we are still able to solve non trivial instances, and find sub-optimal but optimized solutions for the most complex instances. Section IV explains how to encode a sequential pseudo code program description and an architecture description to a set of constraints that an SMT solver can handle. Because the result is a simple Gantt timing diagram, we also detail how to transform back this diagram into an executable binary given the technical limitations of the shared memory heterogeneous platform. Section V presents the obtained experimental benefits, the limits of this approach and how to further stretch this limit with better modeling. We use the running example of a tiled Cholesky algorithm (3x3 to 5x5 tiles) that is built on top of Basic Linear Algebra Subprograms (BLAS). We characterize those Subprograms on each available type of computing resource both in terms of performance and predictability.

II. RELATED WORK

Automatic and assisted parallelization as well as energy optimization are very wide research areas. With the advent of complex computing platforms, carefully modeling the hardware is now a mandatory step. Allocating and scheduling a set of jobs on a set of resources can take different forms in the literature, from very theoretical to very applied. However the three steps of 'characterization, solving, and implementation' can often be identified.

The current practice in consumer electronic devices in order to save energy is to give the user the ability to setup a "policy" or "governor". The default policy is often to set every core to the maximum frequency when a new process is submitted to the device (which is the well known "on-demand" power governor). When the program stops all the cores go back to a low performance (and low energy) state. Because consumer electronic devices must support dynamic task appearance, this is a good approach for a best-effort behavior regarding energy optimization. Embedded signal processing on the contrary often involves a limited set of tasks that must run within a given time budget (such as in a radar processing application), hence the need for static scheduling approaches.

Because every engineer has long been manipulating sequential code in various languages, it is now a difficult problem to extract potential parallelism from this code, and thus previous research has extensively investigated automatic parallelization [START_REF] Pop | Graphite: Polyhedral analyses and optimizations for gcc[END_REF], [START_REF] Bastoul | Putting polyhedral loop transformations to work[END_REF]. Although we give a very naive example of loop nest to graph transformation in order to show the parallelism in the example application, we believe that the reasoning on models in which potential parallelism is less obfuscated (such as graphs and formal models) is also an interesting approach. Pseudo-code or code annotations is also often used for this analysis [START_REF] Bosilca | A generic distributed dag engine for high performance computing[END_REF] and can be interpreted as a first step towards formal models.

Reasoning on formal descriptions of the application (and the architecture) allows to guarantee properties on a system, and is an important trend of the current research. The purpose of those formal models is often to ensure properties of the final design such as safety (absence of deadlocks). Because using the model of a Directed Acyclic Graph of task and messages allows to generate the topological sorts (all the valid sequences), the transformations from the studied formal models to this description is often necessary. Synchronous Dataflow graphs [START_REF] Lee | Synchronous data flow[END_REF] and all its variants give a concise description of parallel applications which allow to validate for safety (absence of deadlocks) and is sometimes used as an input model for those scheduling problems [START_REF] Holmbacka | Energy efficiency and performance management of parallel dataflow applications[END_REF]. Modeling hardware features (such as synchronization or message passing) in the same formal description than the application allows to obtain guarantees on the whole system [START_REF] Millo | Modeling and analyzing dataflow applications on noc-based many-core architectures[END_REF], and not only on the application. In recent work, timed events are also considered in order to check throughput properties of cyclic application graphs mapped on multiple resources [START_REF] Stuijk | Predictable mapping of streaming applications on multiprocessors[END_REF], [START_REF] Stuijk | Sdfˆ3: Sdf for free[END_REF]. Other scheduling techniques include identifying specific task patterns such as fork join, sets of periodic tasks, or pipelines and apply specific heuristics to those patterns [START_REF] Sorel | Syndex: System-level cad software for optimizing distributed real-time embedded systems[END_REF]. Those heuristics sometimes target specific hardware resource organizations such as network on chips [START_REF] Thies | Streamit: A language for streaming applications[END_REF], [START_REF] De Dinechin | A distributed run-time environment for the kalray mppa®-256 integrated manycore processor[END_REF].

Finally, another well known set of techniques involves DES (Discrete Event Simulators) and allows to check the result of a scheduling and mapping strategy when studying a scheduling problem. SystemC is an example of DES that has been used that way, including some of its extensions (TLM) that allow to give a coarse description of the hardware and software models when a RTL or functional implementation is not desired. Some tools are built on top of those DES for fast prototyping (Synopsys MCO Platform Architect [START_REF] Fischer | Power modeling and analysis in early design phases[END_REF]). Those DES often come with verification tools, which allow the creation of testbenches to check specific cases, sometimes including randomization of the inputs. Those tools thus allow estimation but no optimization [START_REF] Lebreton | Power modeling in systemc at transaction level, application to a dvfs architecture[END_REF], [START_REF] Dhanwada | A power estimation methodology for systemc transaction level models[END_REF], unless a model checker is used in place of the DES.

Energy modeling is not often the main purpose of automatic and assisted parallelization [START_REF] Bosilca | A generic distributed dag engine for high performance computing[END_REF], [START_REF] Sorel | Syndex: System-level cad software for optimizing distributed real-time embedded systems[END_REF], [START_REF] Pop | Graphite: Polyhedral analyses and optimizations for gcc[END_REF] although there is a renewed attention for low-power designs. It requires to add non linear expressions in the model and thus raises greatly the complexity for model checking, although simulation is still an option [START_REF] Fischer | Power modeling and analysis in early design phases[END_REF]. Moreover, power modeling involves multiple mechanisms such as frequency scaling and power gating that make both the modeling and the implementation very complex, yet necessary to obtain insightful results [START_REF] Holmbacka | Energy efficiency and performance management of parallel dataflow applications[END_REF].

III. THE ALGORITHM ARCHITECTURE ADEQUATION ISSUE

In this section we show the preliminary work of exhibiting the potential parallelism of the input application model with a pseudo-code to DAG transformation. Then, according to the workflow (Figure 1) we encapsulate the results of the characterization experiments into an application model and an architecture model. In this section we describe what we measure, how, and we give evidence that those models are trustworthy. Because a model is always slightly wrong, we quantify the error that they may introduce. Finally, we describe the energy models that allow us to relate the resources speed and their energy consumption.

A. Application tasks and task graphs

Various input models would fit our needs, as long as there exist a transformation to a graph-like model with a finite number of task and messages. Synchronous dataflow graphs, with the well known Homogeneous SDF transformations is also an eligible candidate, as well as periodic task systems. We achieve this transformation from a pseudo code description of the tiled Cholesky algorithm such as described in the literature [START_REF] Bosilca | A generic distributed dag engine for high performance computing[END_REF], [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF].

1) Pseudo-code to task graph model: Given a number of tiles T of matrix A we use the tiled Cholesky algorithm (Figure 1) to compute the graph of dependencies between the BLAS level-3 functions (GEMM, TRSM, POTRF and SYRK). We obtain the graphs in Figure 2.

Algorithm 1 Pseudo-code for the tiled Cholesky algorithm

Require: number of tiles: T Require: tiled matrix: A for k=0...T do

A[k][k] = POTRF(A[k][k]) for m=k+1...T do A[m][k] = TRSM(A[k][k], A[m][k]) end for for n=k+1...T do A[i][i] = SYRK(A[n][k], A[n][n]) for m=n+1...T do A[n][m] = GEMM(A[m][k],A[n][k],A[m][n]
) end for end for end for Figure 2 shows the data-dependency graph of the tiled Cholesky algorithm for a matrix of 3x3 and 4x4 tiles. We identify separately each function call depending on which tile it reads or writes, hence the numbering in the resulting graph. We add a dependency to the graph for each corresponding read and write. In this model we consider that a task can start only when all its inputs are available.

B. Architectures, Operating Systems and variability

Multi-programming introduces a lot of variability in the execution of data-independent jobs, due for instance to preemption and interrupts. Using a standard environment also comes with some benefits and tools that we can use to check how much variability multi-programming introduces, and even try to avoid it. For instance we expect that thread pinning (force a thread to run and stay on a specific core) will reduce this variability. In this section we describe why variability can be introduced in this setup, and how to keep it low.

WCET (Worst Case Estimation Time) is a difficult topic that we don't fully tackle here. We address the case where the tasks costs do not have much variability which we carefully check in section III-C. We consider the case of an embedded system, where the resources are not shared to multiple users. In that case the variability is mainly due to different behaviors of the tasks depending on their input data. The most flagrant example of "data-dependent" tasks are the sparse matrix operation algorithms, where their dense counterpart are however "dataindependent". For instance we expect a matrix multiplication routine to execute the same number of instructions for fixed matrix sizes whatever values contained in those matrices. This hypothesis does not hold in the case of a sparse matrix multiplication algorithm, because then the algorithm behaves differently depending on the sparsity of the matrices. In that case we could also over-approximate the cost of the matrix multiplication task, which is why the WCET metric is sometimes used. An important part of the characterization step consists in determining the Average Case Estimation Time (ACET), and verifying that its standard deviation has a reasonable value. This reasoning does not apply in the case of an already parallel implementation of a matrix multiplication algorithm, or would require some implementation effort in order to be predictable (which is also what we try to achieve with this methodology).

Our experimental setup allows us to control thread pinning and core frequencies. Because the scheduler of the operating system may decide to switch context or preempt our tasks, we also use a high scheduling priority to further reduce the variability of the task costs. We used the OpenBLAS kernels (implementations of the BLAS level 3 kernels) and carefully check that they run as single threaded kernels, since we aim at building a multithreaded application on top of them. As expected, their multithreaded implementations are much less predictable, with according to our measurements one order of magnitude difference in variability of number of cycles for the matrix multiplication multithreaded implementation without task pinning. This approach distinguishes every sub-program as a sequential execution such that the application that uses it can be both parallel and predictable. A distributed memory architecture such as scratch-pad memories would certainly allow much more predictability, since the implementation would be allowed to control the message passing as well, but because it is not a widespread architecture and because we need also energy saving mechanisms we study this problem on a shared memory architecture, and we try to limit the memory conflicts. The example architecture is a big.LITTLE ARM implementation with frequency scaling illustrated Figure 3. The experiments of sections III-C are conducted with the same experimental setup than the one described in section V-A.

A15 L1 L2 RAM A7 L1 f A7 f A15 A15 L1 L2 RAM A7 L1 f A7 f A15 Fig. 3: bigLITTLE ARM architecture

C. Variability of ground level tasks

We estimate the cost of the task with their number of cycles when executed on the target Odroid-XU3 system. In this section we examine and challenge the low-variability hypothesis. Figure 4 shows the number of cycles when a single threaded BLAS (from the openBLAS benchmark) is executed on the target, along with the Violin error plot and the extreme values. We observe a very stable execution for each of the BLAS we use. We use the coefficient of variation as the variability metric (standard deviation over mean), the results are presented Table I and Figure 4. We perform those measurements by reading the hardware counters through the "perf_events" Linux kernel interface [START_REF] De Melo | The new linux'perf'tools[END_REF]. We use the costs defined in Table I in the next sections in order to obtain parallel schedules for this heterogeneous architecture. In this experiment, the costs of the BLAS kernels are in the same order of magnitude. When the problem combines very small costs with very large costs, the approach may require either to omit the small tasks or to merge them such that the costs of the jobs remain in the same order of magnitudes. Given the Operating Point (OPP) o = (V, f ) of a resource r, we know from previous work [START_REF] Bose | Modeling and analyzing cpu power and performance: Metrics, methods, and abstractions[END_REF] that we can model its instant power with the following:

P (r, o) = C r • V 2 • f • A (1)
We define the variables of Equation 1 and we give their dimension:

• C r [Farads]: the model capacitance of the CPU (roughly proportional to the number of the switching devices, i.e. the number of transistors) • V [Volts]: the voltage of the current operating point • f [Hertz]: the frequency of the current operating point • A [dimensionless]: An activity factor that models the average number of bit flips that the task will cause on the resource. Although it is specific to each pair of task and resource and often approximated to 1. Then for a task t running on the resource r at OPP o, given that the task costs C y cycles when running on this resource, the task will last Cy fo , thus we can model the consumed energy with the expression 2. We observe in this expression that E(t) does not depend on f o , but only on V o (which usually raises with f o ). This means that running a task at a low frequency actually consumes less energy than running it at high frequency, even though it lasts longer. This is the purpose of frequency scaling. Power models sometimes include static power consumptions, which is the amount of power that leaks through the gates as long as the chip is powered. This means that running a task at a very low frequency could actually consume more energy because of static power. Given the measured power levels of Figure 5 we do not consider static power in our approach.

E(t) = C y f o • P (r, o) = C r • C y • V 2 o • A (2)
Figure 5 shows the instant power of the A15 Big resource when it is computing. We used the stress task (from the linux utility "stress", which simply does square root computations) to apply a computation intensive workload. Because for each OPP we know the frequency and the voltage (we set it), and the power sensor allows to sample the instant power, we are able to find realistic values for the capacitance of A7 and A15 cores (resp. 500pF and 100pF). We also sample the instant power when the resource is idle, depending on its running frequency which suggests that the activity factor of the idle state A idle is A active /2. We consider A active = 1 even though this activity factor actually depends on the task which runs and thus might vary depending on the BLAS kernel. We observed small but non negligible variations depending on the task, but the very limited power sensor sampling rate and resolution does not allow to conduct a full characterization. With a proper electronic setting, this value could be integrated into our models with no additional complexity regarding the scheduling problem. 2) Power gating: We consider the static power of the resources but we do not consider that the resources are able to turn on and off during a run. If at least one task has to execute on one resource, then this resource is always on in the steady state system. Then it will be idle in its spare time, at its lowest power OPP but it will still consume power. Thus it is interesting to know if a scheduling problem can be solved by powering ON only a subset of the resources in a machine.

The instant power of the system depends on time because the resources are allowed to change their operating point in time. Given the instant power P (R), the energy with a dimension of W atts • seconds = Joules is thus calculated for each resource r over the duration of the schedule (which means from zero to makespan M ).

E a (r) = t=M t=0 P (r) dt (3)
= C r • t=M t=0 V 2 • f • A dt (4)
Let t o be the opp for task t, this means that each resources r running at opp o consumes:

E(r, o) = C r • t∈T asks|to=o • C y f • V 2 o • f E(r, o) = C r • t∈T asks|to=o C y • V 2 o
Let O r be the set of OPPs of resource r, the energy consumed by the tasks running on r is:

E a (r) = o∈Or E(r, o)
We call that quantity "active" energy consumption.

Finally let P idle be the idle instant power of the resource, t pt the physical time duration of task t and t map the resource where t is mapped. The idle energy can also be expressed as a discrete sum:

E i (r) = P idle • (M - t∈T asks|tmap=r t pt ) (5)
Thus, consumed energy of the system is:

E = sum(E i (r) + E a (r) | r ∈ R) (6)
That E is what we try to minimize. Figure 13 gives a simple example where the consumed energy can be optimized without performance penalty.

3) Motivating example: Lowering the frequency of a task does not always mean that the optimization criteria is degraded. For instance if we want to optimize makespan, we can find very simple application descriptions of tasks and dependencies which give evidences that energy can be optimized with a preserved makespan. Figure 6 is a typical problem instance with a split/join of unbalanced tasks, where slack time can be used to save energy. Generally speaking, this applies to each non critical path. This also gives evidence that the solver will be able to optimize energy if we allow more slack time from the optimal makespan schedule. When the application graph is not a simple unbalanced two workers split/join, it is (computationally speaking) much harder to decide the allocation and the operating point in order to obtain an optimized schedule. For a few tasks with precedences it might even be hard for a human to optimize for makespan, considering only high frequency operating points. Then, the more slack time, the higher possible energy optimization, and the larger state space to explore. Thereafter we explain how to encode those problems into a solver in order to address them automatically. Then we use this constraint programming encoding to schedule the tiles of the tiled matrix Cholesky factorization running on an heterogeneous multicore.

IV. ENCODING INTO CONSTRAINTS OVER FINITE DOMAINS

The metrics that we use to characterize the tasks and the platform actually belong to discrete sets even when some mathematical property is used to define a new variable. For instance the time that a task will last depends on the resource that runs it, and the frequency at which the task runs. This means that the task duration variable is within the set of resources × OPPs. The problem is thus defined as a Constraint Logic Programming over finite domains CLP(fd).

A. Linear Programming

Since the problem mostly includes inequalities and linear equalities such as t stop = t start + t duration , it is a natural idea to model it using linear programming. In heterogeneous scheduling problems, t duration depends on t map , and if multiple operating points exist, it also depends on t opp : we need to introduce a zero-one "decision variable" in both cases. It is easy to encode such constraints using an if-then-else construct, which does not exist in LP most solvers. Another example is the "exclusive OR" constraints of non-overlapping two tasks that use the same resource (see Section IV-C). Some tricks exist but require to introduce artificial variables and quantities that must be chosen carefully. Recent solvers (such as CPLEX) introduce specific constructs that translate to other theories (such as Mixed Integer Programming). Because combining theories is the purpose of "SAT Modulo Theories" solvers, we try to use an SMT solver instead.

B. Costs to finite sets

Let us define the simplest version of the heterogeneous scheduling problem, which is a set of tasks T and a set of resources R. Each resource is associated with a microarchitecture µ ∈ A. The cost of a task in number of cycles is associated to a microarchitecture Cy(t, µ) which is an input of the problem. This means that we can enumerate the possible task costs by iterating the combinations in T × A. Given the obtained Average Case Execution Time estimations of the costs and given the admissible frequencies, the elements of those new set are rationals. Because not every solver have support for the rational type (for instance the SMTLib language [START_REF] Barrett | The smt-lib standard version 2[END_REF] does not define the rational type), we explain how to reduce the set to a set of integers.

For example, given the objective of makespan optimization for the three tasks A,B,C of cost 4MCy, 2MCy, 3MCy on two homogeneous cores (C0 and C1) with different speed 1GHz and 2Ghz: we can bound the problem in time by evaluating one possible schedule, for instance the slowest schedule of 9ms. Then, we know that A is either 4 or 2ms, B: 1 or 2, C: 1.5 or 3. Because we are interested in the solutions where the tasks starts and stops are "aligned", a 1/2 ms which is the Greatest Common Divisor looks like as a very suitable time step. We generalize this by collecting the possible durations of timed events, finding their common divisor and using it as a time step. By construction, this allows to transform every physical time duration into an integer, as well as the time bound. This allows to transform our problem to an Integer Programming problem. Because there exists a finite set of admissible instant powers, the same is conducted for the energy model.

C. Constraints

In this section we list the constraints that we use in order to generate the schedules. Let T be the set of tasks, let R be the set of resources and O be the set of operating points, t map ∈ R and t opp ∈ O. Let M be the set of messages (or precedences), with sender and receiver two functions from M to T . The integer table costs is calculated according to the previous sections. A timebudget is specified, and the makespan is a variable of the problem that may be optimized. Let r ip be the idle power of the resource r, and o p the instant power of opp o. We first define application and allocation scheduling constraints (Algorithm 2), and if the problem requires the power model we add the constraints defined in Algorithm 3.

Algorithm 2 Constraints

Require:

costs table: T × R × O Require: timebudget for t ∈ T do t start , t stop , t duration ∈ N t stop = t start + t duration t stop ≤ makespan makespan ≤ timebudget t duration = costs[t][t map ][t opp ] end for for (v, u) ∈ T 2 | v! = u do u map == v map =⇒ u start ≥ v stop ⊕ u stop ≥ v

D. Back from model to code

Our optimized versions of the tiled Cholesky algorithm use many tools and libraries in order to better control the platform. For the BLAS kernel implementations, the Gnu Scientific Library with OpenBLAS linking was used. libcpufreq allows to change the frequency of a core, libpthread allows to fork and join threads, as well as pinning them to specific cores with the functions defined in sched.h. The result of the SMT solver allows us to build a schedule for each resource, including the data dependencies across resources that we transform to synchronization (because the memory is shared). The functional aspect of the code is tested with assertions that can be removed by the preprocessor with proper symbol definition.

1) Minimal energy schedules: Since the A7 little core can run at the lowest instant power available on the platform, the schedule that minimizes energy is the one that uses only the A7 little core, whatever the resulting makespan. For instance for a 3x3 tiled Cholesky algorithm with minimized energy consumption we obtain the schedule of Figure 7 which is simply one of the many valid topological ordering of the graph presented in Figure 2 POTRF(0,0)->TRSM(1,0)->SYRK(1,0)-> TRSM(2,0)->POTRF(1,1)->GEMM(0,1,2)-> TRSM(2,1)->SYRK(2,0)->SYRK(2,1)->POTRF(2,2) Fig. 7: Topological order of the 3x3 tiled Cholesky task graph 2) Minimal makespan schedules: We now study the schedule obtained when minimizing makespan. Its transformation to a shared memory multithreading environment allows to transform the messages into synchronizations. The highest frequencies (hence the fastest but also the highest instant powers) are selected in order to minimize the makespan, but there is no straightforward solution. The solver outputs a solution (schedule in Figure 8) from which we generate the code skeleton to implement it in the shared memory environment. The messages (arrows in Figure 9) are transformed to POSIX semaphores. Fig. 9: 3x3 Makespan optimal multithreaded pseudo code makespan may contain both synchronizations and frequency changes. We implement the frequency changes with the cpufreq library. This library is the one that implements the energy governors in the linux kernel, and allows to set the governor for each core (userspace, powersave, conservative, ondemand, performance, ...). The typical setup uses either powersave/ondemand/performance. The userspace behavior as its name suggests allows to change the current frequency from a userspace program.

V. RESULTS

In this last section we give the obtained schedules and objective values for different optimization requests and we show how much energy savings one could expect when giving the makespan a bit of slack compared to its optimal value obtained in the previous experiments. Finally, we translate the schedules to executable parallel implementations. Then, we run the obtained binaries on the Exynos5422 platform and we try to evaluate how faithful is the generated implementation.

A. Experimental setup

We run all the experiments on an embedded development board (Odroid XU3) featuring a Samsung Exynos5422 ARM bigLITTLE processor of four A7 cores and four A15 cores such as the one described in Figure 3. The bigLITTLE technology is an heterogeneous MP-SoC architecture designed by ARM. Samsung implemented it in the Exynos5422 Chip, which was used by the Corean company "Hardkernel" to build the Odroid board. The platform runs a Linux Ubuntu 15.10 operating system with a 3.10 kernel that was patched in order to enable hardware cycle counters. In the experiments described in section V we sometimes consider 1,2 or 4 cores of each (big and LITTLE). The platform is equipped with instant power sensors that allow us to sample instant power of each cluster (of big or LITTLE cores) and a performance monitoring unit (PMU) that allows to access relevant metrics such as cycle counts through the perf Linux API [START_REF] De Melo | The new linux'perf'tools[END_REF].

B. Makespan optimization

As expected when the number of tiles grows, the number of tasks and task dependencies also grows. The overall number of samples is constant, thus the number of samples per tile decreases when the number of tiles increases. The costs obtained in Table I have been updated with the corresponding values for the tiles of different sizes corresponding to our experiments. We were able to solve the non trivial examples of 3x3 and 4x4 tiled Cholesky algorithm running on two heterogeneous cores optimally, and we got a sub-optimal (but still optimized) result for the 5x5 tiled Cholesky algorithm. Because the solver is very efficient at finding satisfiable solutions when there is a lot of slack time [START_REF] Gorcitz | On the scalability of constraint solving for static/off-line real-time scheduling[END_REF], even with a much larger number of tasks and messages (resp. 35 tasks and 60 messages in the 5x5 example), it is able to find a solution, and to get an optimal or optimized solution. We experiment different problems with higher theoretical complexity and we report their practical complexity in section V-D.

Because the solver approach also allows to implement heuristics by doing successive calls to the solver and changing the objective function, we can compare the obtained schedules to one of the greedy "As Soon As Possible" (ASAP) schedules. In order to obtain this schedule we sort the tasks by topological order and we add them in this order to the set of constraints, asking for makespan minimization: the strategy becomes greedy, hence sub-optimal but allows as expected to obtain schedules for large instances. We also compare our results to the naive sequential implementation with all tasks running on the fastest resource, at the highest frequency. In either case we obtained significant improvement. For the 5x5 tiled Cholesky algorithm graph, we obtained an optimized but sub-optimal result, still better than naive and ASAP implementations. The results also suggest that the divide and conquer strategy actually works, since the makespan decreases when the number of tiles increases (for the same matrix size). The exhaustive results are presented in Table II.

The synthesized implementation is compared to those predicted values and its predictability is evaluated with the same coefficient of variation used in section III. The results are presented in Table III. Of course running those programs requires a high level of privileges since it requires access to critical hardware features, such as modifying frequencies. This is acceptable on an embedded system since there is probably only one user (sometimes there is even no operating system) but it is probably not acceptable on a shared server with no proper isolation. When the tiles are too large, input data of the BLAS kernels do not fit in the data caches. This can increase variability and lower performance, as the result for the 3x3 case suggests. Adapting tile sizes to cache size such that the generated schedule restricts as much as possible the access to shared levels of the memory would be an interesting addition. The 4x4 and 5x5 results do not suffer from this issue because they compute on smaller tiles. The results in those latter cases are very close to the predicted results, which suggests a very faithful model.

C. Energy optimization

As well as for makespan optimization, the obtained makespans when optimizing for energy are very close to the predicted values (for instance the implementation of Figure 11 lasts 76ms). We now generated Pareto fronts for those models and examine the corresponding implementations. Figure 10 shows a subset of the schedules that are on the Pareto front for the 3x3 algorithm. The power sensor that was used in the section III-D in an on-board device that is connected to the big.LITTLE through an I 2 C Bus. It is not very precise and allows sampling at a slow rate (not at the 1kHz rate that we would need). However, we can at least check the instant power when running different implementations, and the very precise timings allow to calculate the consumed energy. The fastest theoretical schedules for 3 and 4 tiles report a 52.5mJ and 45.1mJ energy consumption. The instant power is sampled when the implementations run and we obtain 1.24W and 1.27W, which corresponds (according to the obtained execution times of Table III) to energy consumption of 66mJ and 51mJ. We do the same for the schedules where the time budget is twice the optimal makespan (Figure 11), those schedules report (resp.) and 30.0 and 30.5mJ, and the implementations report 45.4 mJ and 46.3 mJ. Although the model is not correct, we obtained valid implementations within a given time budget, with a significant energy optimization.

D. Further scalability results

There are a number of direct extensions one can apply to the previous modeling and encoding effort, with in mind the idea of checking how the method can be pushed to its limits in terms of complexity regarding solver efficiency. Already the parametric aspect of Cholesky algorithm allows to scale up a range of applications. Modeling the multiplicity of process frequencies, or multiple cores in each compute block also increases the demand on solvers. We experiment with such extensions now, looking for limitations as they occur. As the Figure 3 suggests, our experimental setup allows to run up to 4 big and 4 little cores. In the previous experiments we have been using only 1+1 cores. We synthesize the results related to scalability in Table III. From the numbers of Table III it appears that, although the scheduling complexity is on the whole related to the size of the problem (the system of constraints to solve), this connection is not as regular as it may seem, and the impact of whether regular or irregular solutions emerge can dominate the size complexity. So certainly much more remains to be done on predicting when (and how) relative dimensioning may favor regularity, and so solver efficiency, in the universe of solutions. Moreover, although in the makespan optimization problem the highest frequencies alone can achieve the optimal results, some of the optimal results can include low frequency tasks (as figure 13 suggests). We observed for instance an impact of the number of operating points on the solver time in the "4 tiles running on 2 little and 2 big cores" problem. 6x6 tiled problems timed out after one hour. The irregularities in the scalability suggest that the constraint programming models can be further improved.

E. Improved modeling for scalability

In the case of the Cholesky matrix factorization, the task graph is not a simple structure. Other applications include more structured patterns such as "split/join" graphs. For instance radars perform multidimensional Fast Fourier Transforms to identify shapes. We experimented with the 2D radar application described in Figure 12. The former approach would consist in instantiating one FFT occurrence for each worker of both "split/join". Another approach consists in instantiating one tasks for each resource, and allow them to have a variable duration as long as their sum complete the total work. The whole problem must also include relevant constraints from Algorithm 2. In that case of an heterogeneous architecture, the latter approach also allows a better load balancing. In the case of four resources (two big and two little cores), the solver returns the optimal solution for up to 6 tasks (n = 6 on Figure 12)for each split/join in a matter of seconds. More complex instances of this naive mode (n > 6) timed out at 10 minutes of search. In this case the naive approach does not scale enough to operate on an image (where usually n > 2 10 ). By defining only two tasks of variable duration, we are able to balance the work of the split/join patterns. Fig. 13: Results for the former (4 tasks for each split join) and latter approach on the radar application VI. CONCLUSION Because shared memory architectures do not allow explicit message passing, the cost of moving data from memory to the L1 caches is not explicitely studied here, it is included in the characterization step of Section III. The constraint programming approach provides a composable environment, thus including message costs for a distributed memory accelerator would be a feasible and interesting addition. It would also be easier to characterize and control the variability in such case. Unfortunately, distributed memory embedded systems such as scratchpad architectures are not widespread.

We do not address the temperature issue in this work. Adding a temperature model would probably require reasoning over a much larger timescale. However, because power consumption raises with temperature this might be an interesting and tricky addition. Future work also includes other model refinement such as introducing other timed events (or penalties) in order to have better modeling which might be necessary for smaller datasets. It is worth noting that refining the models does not always mean that the state space will be larger, for instance introducing a penalty for frequency change actually reduces the state space.

Static scheduling is applicable to other signal or image processing application that admit a divide and conquer algorithm, assuming the ground level tasks are not data-dependent. It is for instance applicable to the Fast Fourier Transform which can be divided into smaller Fast Fourier Transforms, as well as LU and QR matrix factorization and generally speaking any application built on top of the BLAS kernels. In this paper we show with a realistic running example how to use a SMT solver for static scheduling, how to model energy aware architectures and how to trade time for energy with this modeling technique. We investigated the scalability of this approach and applied the same workflow to other examples, including synthetic and realistic applications (FFT, Platooning application). Our results also show that solving realistic structured problems may require more modeling effort. Producing optimized and energy efficient code is a great challenge for the new upcoming battery powered devices.

Fig. 1 :

 1 Fig. 1: Workflow of our approach

Fig. 2 :

 2 Fig.2: Data-dependencies of the tiled Cholesky algorithm for 3x3 and 4x4 tiled matrices (the numbers identify the tile)

Fig. 4 :

 4 Fig. 4: Predictability of BLAS lvl3 kernels on A7 and A15 architectures (in number of cycles), for 200x200 samples tiles

Fig. 5 :

 5 Fig. 5: Instant powers of the tasks vs frequency

Fig. 6 :

 6 Fig.6: Non energy optimized vs optimized split join pattern (red tasks run with high instant power)

Fig. 8 :

 8 Fig.8: Makespan optimal schedule for 3x3 tiled Cholesky algorithm on a 1+1 bigLITTLE architecture (red means high instant power)3) Multiobjective optimization schedules: Finally, the schedules obtained when minimizing energy and bounding

Fig. 10 :

 10 Fig. 10: Pareto front for the 3x3 algorithm

Fig. 11 :

 11 Fig. 11: 4x4 tiled Cholesky schedule optimized (sub-optimal) for energy, with a time budget of twice the optimal

Fig. 12 :

 12 Fig. 12: 2D radar application

TABLE I :

 I Costs of the BLAS kernels (10 6 cycles)

			GEMM SYRK TRSM POTRF
	A15	mean c v (%)	11 0.29	6 0.21	7 0.19	9 0.08
	A7	mean c v (%)	39 0.11	20 0.23	21 0.11	16 0.64
	we found realistic values for C r and activity factor on the
	big.LITTLE Exynos5422 target.
	1) Dynamic Voltage and Frequency Scaling:

  Energy related variables and constraints for o ∈ O do o uptime = t∈T (t duration | t map == r, t opp == o) o energy = o p * o uptime end for for r ∈ R do r uptime = o∈O (o uptime ) r downtime = makespan -r uptime if r uptime == 0 then r energy = 0 else r energy = o∈O (o energy ) + r ip × r downtime end if end for

	start
	end for
	for m ∈ M do
	receiver(m) start ≥ sender(m) stop
	end for
	Algorithm 3

TABLE II :

 II Results obtained with the modeling and solving (time in milliseconds)

		3x3	4x4	5x5
	N.of tasks (and messages) 10 (12) 20 (30) 35 (60)
	Optimal makespan	46	39	≤ 38
	Sequential makespan	55	51	49
	ASAP makespan	55	45	44
	Energy optimal makespan 350	349	343

TABLE III :

 III Execution time and variability of the generated implementations

		3x3 4x4 5x5
	Time (ms) 53.2 40.5 39.6
	c v (%)	0.89 0.53 0.77

TABLE IV :

 IV Solver time vs number of cores and number of tiles

			#cores (big+LITTLE)
			1+1	2+2 4+4
		3x3	2s	2s	3s
	#tiles	4x4 1min48s	8s	11s
		5x5	>1h	10s 28s


D Galassi
email: davide.galassi@univ-amu.fr

P Tamain

C Baudoin

H Bufferand

G Ciraolo

N Fedorczak

Ph Ghendrih

N Nace

E Serre

Flux expansion effect on turbulent transport in 3D global simulations

The flux expansion effect on the Scrape-Off Layer equilibrium is inspected through TOKAM3X 3D turbulence simulations. Three magnetic equilibria with analytically controlled flux expansion are built, representing respectively a positive, a null and a negative Shafranov shift. Turbulent E × B fluxes across flux surfaces show similar amplitudes and poloidal distributions in all cases. The ballooning nature of the interchange instability is recovered, with an enhancement of turbulence in the vicinity of the limiter, probably due to a Kelvin-Helmoltz instability. Interestingly, the flux surface averaged density decay length is found to be shorter almost by a factor 2 in the case of flux surfaces compressed at the Low Field Side midplane, indicating the presence of unfavourable conditions for the turbulent transport. The difference in the magnetic field line shape is pointed out as a mechanism which affects the turbulent transport across the flux surfaces. Indeed the unstable region has a larger parallel extension when flux expansion in the Low Field Side is larger. Moreover, the configuration with a lower magnetic shear at the Low Field Side (LFS) midplane 1 shows a more unstable behaviour. The role of this parameter in turbulence stabilization is qualitatively evaluated. The difference in the distribution of transport along the parallel direction is shown to affect also the parallel flows, which are analyzed for the three proposed cases.

Introduction

Recently, the fluid turbulence code TOKAM3X has been run in a divertor configuration, in order to analyze the effects of turbulent transport on global flows in a X-point geometry [START_REF] Galassi | Drive of parallel flows by turbulence and large-scale e x b transverse transport in divertor geometry[END_REF]. One peculiarity of the turbulent transport arising from this simulation, is the fact that the turbulent fluxes tend to adapt to the tokamak magnetic geometry. It is shown, indeed, that the cross-field turbulent flux is enhanced where the flux expansion is stronger, so at the top of the machine and at the X-point. Turbulent structures indeed, at a first order, tend to be damped when they are not field aligned. In order to follow the flux surface expansion so, turbulent structures must propagate with different velocities, according to the flux tube radial extension in the physical space. This distance is regulated by the flux expansion, which has thus a direct geometric effect on turbulent fluxes. Putting ourselves in a framework independent from the magnetic flux expansion, we can actually study the locations where the turbulent transport is more effective. Simulations in divertor configurations, for example, have shown that the LFS midplane is actually the location where the turbulent transport across the flux surfaces is more effective. The scope of this paper is to inspect more deeply the effect of flux expansion on turbulent transport, independently form the direct geometrical effect on cross-field fluxes just described. In other words, we aim at studying the effect of a geometric variation in the real space on the turbulence studied in the magnetic space. This work is carried out through numerical simulations with the 3D fluid turbulence code TOKAM3X [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF], whose flexibility in the geometry definition allows to model a large variety of magnetic configurations. Some works on this subject have been carried out with other codes considering a single flux tube geometry and averaging the dynamics of seeded blobs [START_REF] Walkden | Characterization of 3D filament dynamics in a MAST SOL flux tube geometry[END_REF]. A global approach is instead retained here, with turbulence arising naturally in the plasma edge, caused mainly by the interchange instability, and in the framework of a realistic three-dimensional geometry.

The TOKAM3X code

The set of equations solved by TOKAM3X, and the associated boundary conditions, can be found in [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF]. It is a drift-reduced, electrostatic turbulence model. For the following analysis, TOKAM3X is run in isothermal mode, with a constant temperature for both species, T e = T i = T 0 , where T 0 is the reference temperature. This means that in our case the plasma density coincides with the static pressure. The electric drift velocity is expressed as u E = ( B × ∇Φ)/B 2 , while the curvature drift velocity is expressed as u i/e ∇B = ±2T i/e ( B × ∇B)/B 3 , with the sign + for the ions i and -for the electrons e. Neutral dynamics is not taken into account. This is an important limitation in the study of the SOL physics. Nevertheless, we are trying to understand a basic physical mechanism as the turbulent transport, so not including the neutral dynamics could help us studying this problem at a basic level, being independent from other physical mechanisms. This approach, moreover, can be representative in situations where the neutral dynamics has a lower impact, as for example in the sheath-limited regime. Four parameters intervene in the physical model. Three of them are diffusion coefficients, for density, parallel momentum and vorticity, and are set for the simulation as D N,Γ,W = 0.5 • 10 -2 ρ 2

L ω c , where ρ L is the ion Larmor radius and ω c the ion gyrofrequency. These coefficients determine the efficiency of the collisional transport, and help damping the turbulent structures having a scale length comparable to the mesh grid. The fourth parameter is the parallel resistivity, and it is set to η = 1 • 10 -5 B 0 /(en 0 ), where n 0 is the reference density used for the normalization. The parallel resistivity determines the losses in the parallel direction, thus affecting the SOL operational regime. The particle source has a Gaussian shape in the radial direction, with a halfwidth of a/16, where a is the minor radius. It is poloidally and toroidally constant, and it is located at the inner boundary of the domain, in the closed flux surfaces region.

Simulation in limiter configuration with flux surface shift

We want to inspect the role of flux expansion on turbulent transport and on global equilibrium. Since we aim at studying the flux expansion effect at a basic level, our scope is to disentangle it from all the other possible geometrical aspects related to the divetor configuration, as, for example, the poloidal field singularity at the X-point and the presence of the Private Flux Region. For this reason, three analytic magnetic equilibria are built.

Simulated geometry and conditions

A toroidally limited magnetic configuration is used, where the limiter is modelled as infinitely thin and localized poloidally at the bottom of the machine.

The minor radius a of the simulated tokamak measures 256 ρ L , and the aspect ratio is 3.4. The three magnetic equilibria differ in the flux surface shifts. Flux surfaces are circular, but their centre horizontal position R ψ is shifted according to the following equation:

R ψ (ψ) = R 0 + ∆(ψ) = R 0 + S • (r(ψ) -a) (1)
Where ∆ is the Shafranov shift, ψ is the flux surface coordinate, and the shift parameter S for the three simulations is respectively 1/3, 0 and -1/3. These three cases will be named in the following "inner shift", "no shift"and "outer shift"respectively. We define the flux expansion as:

f x (ψ, θ) = < ∇ψ (ψ, θ) > θ ∇ψ (ψ, θ) (2)
With this definition and the chosen parameters for the shift, in the vicinity of the separatrix the value of the flux expansion on the expanded side, in the shifted cases, is twice the value on the compressed side and 4/3 times the value in the case without shift. The poloidal flux function Ψ is adjusted in order to obtain the same profile of mean security factor q(ψ) for the three simulations. In particular, q(ψ) has a parabolic shape, with a value around 4 at the separatrix. The toroidal magnetic field is chosen to vary as 1/R, with a normalized unitary value at the magnetic axis. As a result, the poloidal field has the shape shown in Figure 1.

One can notice from figure 1 that while in the inner shift case the poloidal field has its maximum at the HFS midplane, for the outer shift case this maximum is located at the LFS midplane. The total connection length is slightly varied in the three configurations. In particular the outer shift geometry has a connection length at the separatrix around 5% shorter than the case without shift, while for the inner shift configuration this value is roughly 5% bigger.

Simulation results

In order to characterize the turbulent transport in the presented geometry, we start by analyzing the turbulent E × B fluxes across the flux surfaces. Making an averaging operation on time and on the toroidal direction:

< N u ψ E > t,ϕ = < Ñ ũψ E > t,ϕ + < N > t,ϕ < u ψ E > t,ϕ (3)
Turbulent fluxes, identified by the second term in the RHS of (3), are associated with the fluctuating components of the density and electric potential fields. In contrast with the simulations in divertor geometry, here the turbulent fluxes are dominant with respect to the mean-field fluxes, the first term of the RHS in (3), associated instead on the average components of the fields. Mean-field fluxes are important only at the vicinity of the targets, where the sheath drives strong potential gradients in the poloidal direction. In order to investigate the turbulence transport efficiency, we analyze the fluxes measured in the magnetic space, thus normalized by the local value of the flux expansion, as defined in (2). We can notice from figure 2 that in the poloidal direction, turbulent flux distribution does not vary appreciably, and it is mainly centred at the LFS midplane. The geometrical effect of flux expansion, already found in divertor geometry, is confirmed: if the velocities were measured in the physical space, we would find values proportional to the local flux expansion. In the inner shift case so, fluxes at LFS midplane in the ψ direction would be, at a first order, larger by a factor two with respect to the outer shift configuration. Studying the fluxes independently from the flux expansion allow us to understand where the turbulent transport is more efficient in crossing the flux surfaces. Flux analysis in the magnetic space shows indeed the ballooning-like character of the turbulent flux, typical of the interchange instability, which dominates in this region over the other possible instability mechanisms (Kelvin-Helmoltz, drift waves) for the set of parameter used in these simulations. A considerable outwards turbulent flux, can be noticed in the three cases in the vicinity of the limiter on the LFS. This turbulence enhancement is likely to be due to a Kelvin-Helmoltz instability, caused by the parallel shear flow between the SOL and the edge region [START_REF] Schwander | Parallel shear flow instability in the tokamak edge[END_REF]. Indeed, while the parallel flow is almost null in the closed flux surface region, the parallel Mach number passes brutally to nearly sonic values when crossing the separatrix. The instability criterion expressed in [START_REF] Schwander | Parallel shear flow instability in the tokamak edge[END_REF] is fulfilled, and the Kelvin-Helmoltz process seems to co-exist with the interchange at the outer side of the limiter in the Scrape-Off Layer (SOL). Globally, the cross-field turbulent flux does not change appreciably in the three cases. This fact is not astonishing, if we think that our system is fluxdriven, and the particle source is fixed and equal in the three simulations. Interestingly, simulations show that the equilibrium in the radial direction changes in the three cases. Since the cross-field flux is somehow fixed by the system, the density gradient adapts according to the efficiency of the turbulent transport. Figure 3 shows radial profiles for the average density in the three cases. We notice that density gradients in the magnetic space get steeper for bigger flux compression on the LFS. The case without shift is in an intermediate situation, as could be expected, but it seems actually closer to the inner shift case. An exponential fit on the radial profiles in the magnetic space, shown above, gives the following results for the radial decay length:

λ IS N = 56 λ N S N = 47 λ OS N = 34
Where superscripts stand for "Inner Shift", "No Shift"and "Outer Shift". Different possible mechanisms can contribute to the turbulence stabilization.

Let's first examine the role of the collisional transport. The most interesting region to focus on is the LFS midplane, where the majority of the turbulent cross-field flux is concentrated. At this position, for the outer shift case, the flux surfaces are more compressed by a factor two with respect to the inner shift case. In our physical model, the destabilizing terms are linked to the E × B and the diamagnetic advections, which appear under a divergence form. The interchange instability is thus driven by terms under an order one derivative. The collisional transport is modelled with a laplacian operator and the diffusion coefficients determine its weight with respect to the other terms. We use in these simulations a spatially homogeneous diffusion coefficient. Thus when compressing the flux surfaces, looking at the physical space, the destabilizing terms grow proportionally, while the diffusive ones grow with a quadratic trend. Thus the effectiveness of the collisions grows with the flux compression, and a steeper gradient in the magnetic space could be needed to drive the turbulence. In order to evaluate this effect, two simulations have been run with a simple limiter geometry, varying the diffusion coefficient. This procedure has been chosen in order to simulate artificially the local increase (or decrease, depending on the flux expansion value) of diffusion coefficient described above. In the first simulation the diffusion coefficient has been divided by a factor 5, while in the second it has been multiplied by a factor 5. This modulation of the diffusion coefficient value has not put in evidence any remarkable changes in the average density profiles. This result is reassuring, since one of the hypothesis underneath the TOKAM3X physical model consists in the negligibility of the diffusive transport with respect to the convective one. We are now sure that the simulations with the three different magnetic equilibria are in a parameter region where the collisional transport has a negligible impact. However, the interpolation scheme used in the code could induce a small numerical diffusion whose evaluation is more challenging, and cannot be totally excluded. The second mechanism involves the shape of the field lines in the three different cases. From a linear analysis done on the physical system of equations solved by the code [START_REF] Tamain | Etude des flux de matière dans le plasma de bord des tokamaks : Alimentation, transport et turbulence[END_REF], we know that the destabilizing term for the interchange instability is:

g = 4 ∂ ∂ r 1 B ( 4
)
The system is unstable when g is positive, so when the local curvature is favourable for the interchange mechanism. If we want to explain the trends of the average density decay length, we need however to calculate an integrated value which characterizes the average turbulent transport over a flux surface. For this reason, we calculate the integral over a poloidal turn in the parallel direction:

< g > θ = 2π 0 g B B p dθ / 2π 0 B B p dθ ( 5
)
Calculating the value of < g > θ in the three cases gives the following results, for the first flux outside the separatrix, we obtain:

< g > IS θ = 1.6 • 10 -3 < g > N S θ = 9.6 • 10 -5 < g > OS θ = -1.4 • 10 -3
The calculated values vary along the radial direction in a similar way for the three cases, so that the difference between them is almost constant. These results confirm us that the plasma edge is more unstable in the inner shift configuration. Indeed, the g term, which has a cosine shape in θ direction, is shifted towards positive values when the shift parameter S is positive. In our case this happens for the inner shift equilibrium. Here the region with a positive g term is the most extended in the poloidal direction, as a direct effect of the positive shift. Moreover, the parallel length in the unstable region is even longer for the inner shift case, in which the local security factor is higher at the LFS. This means that globally, a particle will cover a longer path along the field line in an unstable region, so it will be more affected by the cross-field turbulent transport. Indeed, in the case of inner shift, smaller density gradients are needed to drive the same turbulent fluxes, so we can say that it is the most unstable. The global g parameter shows marked differences between the three cases, so it could, at least partially, explain the macroscopic differences among them. The different magnetic shear in the three situation is identified as another possible mechanism which affects the turbulent transport. As mentioned before, the profile of the security factor q(ψ) is maintained constant over the three simulations. Nevertheless, looking at the dependency of q on the radial coordinate, this equality is not valid any more, since the spacing between the flux surfaces is not the same. In particular, the q(r) profile is much steeper for the outer shift cases. This changes, by consequence, the local magnetic shear. In a toroidal geometry with a circular cross-section, the magnetic shear can be defined as:

s(r) = r q(r) dq(r) dr (6)
Where q is defined as in [START_REF] Wesson | Tokamaks[END_REF]. Calculating the shear value at the LFS midplane, for a radial position correspondent to the separatrix, we obtain:

s IS = 1.6 s N S = 2.1 s OS = 3.3
Also this parameter so varies strongly among the three simulations. A magnetic shear stabilization is reported experimentally in cases of negative magnetic shear [START_REF] Antonsen | Physical mechanism of enhanced stability from negative shear in tokamaks: Implications for edge transport and the L-H transition[END_REF]. The shearing mechanism due to the magnetic field depends (as explained in [START_REF] Antonsen | Physical mechanism of enhanced stability from negative shear in tokamaks: Implications for edge transport and the L-H transition[END_REF]) on the alignment between the blobs and the horizontal direction, along which the interchange mechanisms acts more efficiently. This is a consequence of the fact that the charge separation happens in the vertical direction, caused by the ∇B drift. The more the density structures stay aligned with the horizontal direction, the more the turbulent transport will be effective. From this simple explanation we can deduce that the shearing efficiency has a minimum for a given magnetic shear value. Both linear and non-linear gyrokinetic simulations in simplified geometry have confirmed this trend [START_REF] Waltz | A gyro-Landau-fluid transport model[END_REF], [START_REF] Bourdelle | ITPA Group, and the International ITB Database Work. Impact of the α parameter on the microstability of internal transport barriers[END_REF] for ITG instability. Analyzing the shape of the blobs near the separatrix in our simulations, we find turbulent structures which get more and more tilted with respect to the horizontal direction for increasing magnetic shear. In particular, in the cases with a more elevated shear, the structures tends to stay aligned to the horizontal direction only for a smaller poloidal span. It is plausible so, that the elevated shear values which characterize our turbulence simulations are larger than the optimum value for the turbulence development, and so that a higher magnetic shear leads to an improved stability. This picture is coherent with our simulation results. The magnetic shear effect, nevertheless, is not limited to a purely geometric feature. Indeed it is known to act in synergy with the E ×B shear, as pointed out in [START_REF] Fedorczak | Dynamics of tilted eddies in a transversal flow at the edge of tokamak plasmas and the consequences for LH transition[END_REF]. This analysis would thus need a further scan in the magnetic shear parameter for the three shown cases, in order to build the stability curve appropriately. This goes beyond the scopes of the present paper.

SOL paralell flow analysis

An analysis of the parallel flows in the SOL is useful to have a more complete sight on the global equilibrium. Figure 4 shows the poloidal profiles of the parallel Mach number resulting from the simulations, at a flux surface in the middle of the SOL. We can notice from figure 4 that the trends in parallel Mach number poloidal profiles are very similar. In particular, the stagnation point is very close to the LFS midplane and almost coincides for the three cases. According to the model presented in [START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF], this tells us that the center of mass of the particle source in the SOL is placed at the LFS midplane. This confirms the fact that the interchange turbulence is the main driver for the parallel flow in the SOL. Mach profile in the outer shift case shows a slight difference with respect to the other cases, mainly far from the stagnation point. In particular, the parallel Mach number is higher in absolute value in these regions. From the model presented in [START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF], we know that what determines the Mach number is the distribution along the parallel direction of the particle source. Since the security factor, as explained before, is lower for the outer shift configuration, a fixed poloidal displacement results in a smaller path in the parallel direction. So even if the poloidal distribution of the cross-field flux, which is the particle source for the parallel flows, is similar, its profile along the parallel direction is different. More specifically, the particle source results more localized in the parallel direction with respect to the other geometries, resulting in higher Mach number values far from the stagnation point.

Conclusion

TOKAM3X 3D global turbulence simulations are run using three magnetic equilibria characterized by different flux expansions at the LFS midplane. Turbulent fluxes across the flux surfaces have been studied, referring to magnetic coordinates which are independent from the flux expansion itself. The resulting fluxes have a similar poloidal distribution, due to the dominant interchange mechanism, and a comparable amplitude, due to the fact that the system is driven by a fixed source. Steeper density gradients along the magnetic coordinate are found in the case of outer shift, so when flux surfaces are more compressed at the LFS midplane, denoting an improvement in turbulence stabilization. The case without shift presents an intermediate situation, and the case where the maximum of the flux expansion is situated at the LFS midplane shows the flattest gradients. In order to explain this, the curvature term g, which locally drives the turbulence, is averaged over the poloidal direction. Results shows that, in the case with flux expansion concentrated at the midplane, particles cover a longer parallel path in a favourable region for the instability: this comes both from the shift of the flux surfaces itself and from the higher security factor at the LFS midplane. Macroscopic differences are also found in the magnetic shear parameter, calculated at the LFS midplane. A stronger magnetic shear in the case of weaker flux expansion at LFS midplane is believed to help in the dissipation of the turbulent structures, improving thus the stabilization. Globally, the simulations disclose the importance of the magnetic geometry in the determination of the turbulent cross-field transport. There is an interesting interplay between features linked to the physical space and to the magnetic one, since the variations in the real space of magnetic features, as the toroidal field or the security factor, play a fundamental role in driving the turbulence. A variation in the flux expansion affects both the geometrical and magnetic characteristics, thus conditioning the overall equilibrium. Results are in line with similar simulations found in literature, that already underlined the importance of flux expansion and magnetic shear for the blob propagation [START_REF] Walkden | Characterization of 3D filament dynamics in a MAST SOL flux tube geometry[END_REF]. In divertor geometry, the described effects are likely to be greatly enhanced. In this magnetic configuration indeed, the flux expansion presents a strong poloidal variability, notably between the LFS midplane and the X-point region. Higher values of magnetic shear at the midplane, as defined in the previous section, are achieved, because of the presence of the X-point. Taking into account all these aspects, the divertor configuration should lead to to a higher turbulence stabilization, taking to extreme the case of the outer shift. The presented work, putting in evidence the effects of the flux expansion on turbulent transport, will allow thus to better understand the simulations run in divertor configuration.

Figure 1 :

 1 Figure 1: Poloidal magnetic field in the three simulated cases: a) inner shift b) no shift c) outer shift.

Figure 2 :

 2 Figure 2: 2D poloidal maps of turbulent flux across the flux surfaces. From left to right: inner shift, no shift and outer shift cases

Figure 3 :

 3 Figure 3: Radial profile of density averaged over the flux surface, remapped on reference flux surfaces.

Figure 4 :

 4 Figure 4: Poloidal profile of the parallel Mach number at a specific flux surface in the Scrape-Off Layer

Acknowledgements

This work has been carried within the framework of EUROfusion Consortium and has received funding from the EURATOM research and training programs under grant agreement WP14/ER/CEA09. The authors would like to thank the financial support by the AMIDEX project KFC. The view and opinions expressed herein do not necessarily reflect those of the European Commission.