Emmanuelle Saillard
email: emmanuelle.saillard.ocre@cea.fr

Hugo Brunie
email: hugo.brunie.ocre@cea.fr

Patrick Carribault
email: patrick.carribault@cea.fr

Denis Barthou
email: de-nis.barthou@labri.fr

PARCOACH Extension for Hybrid Applications with Interprocedural Analysis

Supercomputers are rapidly evolving with now millions of processing units, posing the questions of their programmability. Despite the emergence of more widespread and functional programming models, developing correct and effective parallel applications still remains a complex task. Although debugging solutions have emerged to address this issue, they often come with restrictions. Furthermore, programming model evolutions stress the requirement for a validation tool able to handle hybrid applications. Indeed, as current scientific applications mainly rely on MPI (Message-Passing Interface), new hardwares designed with a larger node-level parallelism advocate for an MPI+X solution with X a shared-memory model like OpenMP. But integrating two different approaches inside the same application can be error-prone leading to complex bugs. In an MPI+X program, not only the correctness of MPI should be ensured but also its interactions with the multi-threaded model. For example, identical MPI collective operations cannot be performed by multiple non-synchronized threads. In this paper, we present an extension of the PARallel COntrol flow Anomaly CHecker (PARCOACH) to enable verification of hybrid HPC applications. Relying on a GCC plugin that combines static and dynamic analysis, the first pass statically verifies the thread level required by an MPI+OpenMP application and outlines execution paths leading to potential deadlocks. Based on this analysis, the code is selectively instrumented, displaying an error and interrupting all processes if the actual scheduling leads to a deadlock situation.

Introduction

The evolution of supercomputers to Exascale systems raises the issue of choosing the right parallel programming models for applications. Currently, most HPC applications are based on MPI. But the hardware evolution of increasing core counts per node leads to a mix of MPI with shared-memory approaches like OpenMP. However merging two parallel programming models within the same application requires full interoperability between these models and makes the debugging task more challenging. Therefore, there is a need for tools able to identify functional bugs as early as possible during the development cycle. To tackle this issue, we designed the PARallel COntrol flow Anomaly CHecker (PARCOACH) that combines static and dynamic analyses to enable an early detection of bugs in parallel applications. With the help of a compiler pass, PARCOACH can extract potential parallel deadlocks related to control-flow divergence and issue warnings during the compilation. Not only the parallel contructs involved in the deadlock are identified and printed during the compilation, but the statements responsible for the control-flow divergence are also outputed. In this paper, we propose an extension of PARCOACH to hybrid MPI+OpenMP applications and an interprocedural analysis to improve the bug detection through a whole program. This work is based on [START_REF] Saillard | MPI Thread-Level Checking for MPI+OpenMP Applications[END_REF] and extends [START_REF] Saillard | Static/Dynamic Validation of MPI Collective Communications in Multi-threaded Context[END_REF] with more details and an interprocedural analysis. To the best of our knowledge, only Marmot [START_REF] Hilbrich | Detection of Violations to the MPI Standard in Hybrid OpenMP/MPI Applications[END_REF] is able to detect errors in MPI+OpenMP programs. But as a dynamic tool, Marmot detects errors during the execution and is limited to the dynamic parallel schedule and only detects errors occuring for a given inputset whereas our approach allows for static bug detection with runtime support and detects bugs for all possible values of inputs.

In the following we assume that all programs are SPMD MPI programs and all MPI collective operations are called with compatible arguments (only the MPI COMM WORLD communicator is supported). Therefore, each MPI task can have a different control flow within functions, but it goes through the same functions for communications. Issues related to MPI arguments can be tested through other tools.

Motivating Examples

The MPI specification requires that all MPI processes call the same collective operations (blocking and non-blocking since MPI-3) in the same order [START_REF]Message Passing Interface Forum[END_REF]. These calls do not have to occur at the same line of source code, but the dynamic sequence of collectives should be the same otherwise a deadlock can occur. In addition, MPI calls should be cautiously located in multi-threaded regions. Focusing only on MPI, in Listing 1, because of the conditional in line 2 (if statement), some processes may call the MPI Reduce function while others may not. Similarly, in Listing 2, some MPI ranks may perform a blocking barrier (MPI Barrier) while others will call a non-blocking one (MPI Ibarrier). The sequence is the same (call to one barrier), but this blocking/non-blocking matching is forbidden by the MPI specification. Regarding hybrid MPI+OpenMP applications, the MPI API defines four levels of thread support to indicate how threads should interact with MPI: MPI THREAD SINGLE, MPI THREAD FUNNELED, MPI THREAD SERIALIZED and MPI THREAD MULT IPLE. MPI processes can be multithreaded but the MPI standard specifies that "it is the user responsibility to prevent races when threads within the same application post conflicting communication calls" [START_REF]Message Passing Interface Forum[END_REF]. In Listing 2, MPI calls are executed outside the multithreaded region. This piece of code is therefore compliant with the MPI THREAD SINGLE level. But MPI communications may appear inside OpenMP blocks. For example, the MPI point-to-point function at line 7 in Listing 3 is inside a master block. The minimum thread level required for this code is therefore MPI THREAD FUNNELED. However, calls located inside a single or master block may lead to different thread support. Indeed, in Listing 4, two MPI Reduce are in different single regions. Because of the nowait clause on the first single region, these calls are performed simultaneously by different threads. This example requires the maximum thread support level i.e., MPI THREAD MULTIPLE.

These simple examples illustrate the difficulty for a developer to ensure that MPI calls are correctly used inside an hybrid MPI+OpenMP application. A tool able to check, for each MPI call, in which thread context it can be performed would help the application developer to know which thread-level an application requires. Furthermore, beyond this support, checking deadlock of MPI collective communications in presence of OpenMP constructs can be very tricky. In this paper, we propose an extension of PARCOACH to tackle these issues, with the help of an interprocedural analysis to improve the compile-time detection.

Section 2 gives an overview of the PARCOACH platform with a description of its static and dynamic analyses for hybrid MPI+OpenMP applications. Then, Section 3 describes an interprocedural extension of the PARCOACH static pass. Section 4 presents experimental results and finally Section 5 concludes the paper.

PARCOACH Static and Dynamic Analyses for Hybrid Applications

PARCOACH uses a two-step method to verify MPI+OpenMP applications as shown in Figure 2. The first analysis is located in the middle of the compilation chain, where the code is represented as an intermediate form. Each function of a program is depicted by a graph representation called Control Flow Graph (CFG). PARCOACH analyses the CFG of each function to detect potential errors or deadlocks in a program. When a potential deadlock is detected, PARCOACH reports a warning with precise information about the possible deadlock (line and name of the guilty MPI communications, and line of conditionals responsible for the deadlock). Then the warnings are confirmed by a static instrumentation of the code. Note that whenever the compile-time analysis is able to statically prove the correctness of a function, no code is inserted in the program, reducing the impact of our transformation on the execution time. If deadlocks are about to occur at runtime, the program is stopped and PARCOACH returns error messages with compilation information. This section describes the following new features of PARCOACH: (i) detection of the minimal MPI thread-level support required by an MPI+OpenMP application (see [START_REF] Saillard | MPI Thread-Level Checking for MPI+OpenMP Applications[END_REF] for more details) and (ii) checking misuse of MPI blocking and nonblocking collectives in a multi-threaded context (extension of [START_REF] Saillard | Static/Dynamic Validation of MPI Collective Communications in Multi-threaded Context[END_REF]).

MPI Thread-Level Checking

This analysis finds the right MPI thread-level support to be used and identifies code fragments that may prevent conformance to a given level. Verifying the compliance of an MPI thread level in MPI+OpenMP code resorts to check the placement of MPI calls. To determine the thread context in which MPI calls are performed, we augment the CFGs by marking the nodes containing MPI calls (point-to-point Fig. 2 PARCOACH two-step analysis overview and collective). Then, with a depth-first search traversal, we associate a parallelism word to each node. As defined in [START_REF] Saillard | MPI Thread-Level Checking for MPI+OpenMP Applications[END_REF], a parallelism word is the sequence of OpenMP parallel constructs (P:parallel, S:single, M:master and B:barrier for implicit and explicit barriers) surrounding a node from the beginning of the function to the node. The analysis detects CFG nodes containing MPI calls associated to parallelism words defining a multithreaded context and forbidden concurrent calls. Based on this analysis, the following section describes how collectives operations can be verified in a multithreaded context.

MPI Collective Communication Verification

This analysis proposes a solution to check the sequence of collective communications inside MPI+OpenMP programs. PARCOACH verifies that there is a total order between the MPI collective calls within each process and it ensures that this order is the same for all MPI ranks. Our analysis relies on checking 3 rules:

1. Within an MPI process, all collectives are executed in a monothreaded context; 2. Within an MPI process, two collective executions are sequentially ordered, either because they belong to the same monothreaded region or because they are separated by a thread synchronization (no concurrent monothreaded regions); 3. The sequence of collectives are the same for all MPI processes (i.e., sequences do not depend on the control flow).

A function is then said to be potentially statically incorrect if at least one of the three categories presented in Figure 3 is verified. This section describes how these error categories can be detected. Category 1 Detection: This phase of the static analysis corresponds to the detection of MPI collectives that are not executed in a monothreaded region. To this end, we use the parallelism words defined in [START_REF] Saillard | Static/Dynamic Validation of MPI Collective Communications in Multi-threaded Context[END_REF]. A parallelism word defines a monothreaded context if it ends with an S or an M (ignoring Bs). If the parallelism word has a sequence of two or more P with no S or M in-between, it implies the parallelism is nested. Even if the word ends with an S or M, one thread for each thread team can execute the MPI collectives. For this part, it is not necessary to seperate single from master regions. So the finite-state automaton in [START_REF] Saillard | MPI Thread-Level Checking for MPI+OpenMP Applications[END_REF] is simplified into the automaton presented Figure 4. It recognizes the language of parallelism words corresponding to monothreaded regions. States 0 and 2 are the accepting states and the language L defined by L = (S|M|B|PB * S|PB * M) * contains the accepted words (parallelism words ending by S or M without a repeated sequence of P).

Category 2 Detection: For this analysis, MPI collective operations are assumed to be called in monothreaded regions, as defined in the previous section. However, different MPI collectives can still be executed simultaneously if monothreaded regions are executed in parallel. This phase corresponds to the detection of MPI collective calls in concurrent monothreaded regions.

Two nodes n 1 and n 2 are said to be in concurrent monothreaded regions if they are in monothreaded regions and if their parallelism words pw[n 1 ] and pw[n 2 ] are respectively equal to wS j u and wS k v where w is a common prefix (possibly empty) with j = k, u and v words in (P|S|B) * .

Category 3 Detection: Once the sequence of MPI collective calls are verified in each MPI process, we must check that all sequences are the same for all processes.

To verify that we rely on Algorithm 1 proposed in [START_REF] Saillard | PARCOACH: Combining static and dynamic validation of MPI collective communications[END_REF] with the extension of nonblocking collectives detailed in [START_REF] Jaeger | Correctness Analysis of MPI-3 Non-Blocking Communications in PARCOACH[END_REF]. It detects MPI blocking and non-blocking collective mismatches by identifying conditionals potentially leading to a deadlock situation (set S). A warning is also issued for collective calls located in a loop as they can be called different times if the number of iterations is not the same for all MPI processes. Then MULTITHREADED REGIONS and CONCURRENT CALLS procedures respectively detect categories 1 and 2 of errors. Finally the STATIC PASS procedure detects category 3 of errors.

Static Instrumentation

The compile-time verification outputs warnings for MPI collective operations that may lead to an error or deadlock. Nevertheless the static analysis could lead to false positives if the actual control-flow divergence is not happening during the execution. To deal with this issue, we present a dynamic instrumentation that verifies warnings emitted at compile-time.

Algorithm 2 Library Functions To Check MPI collectives To dynamically verify the total order of MPI collective sequences in each MPI process, validation functions (CCipw and CCcc) are inserted in nodes in the sets Sipw and Scc generated by the static pass (see Algorithm 1). These functions are depicted in Algorithm 2. Function CCipw detects incorrect execution parallelism words and Function CCcc detects concurrent collective calls. To dynamically verify the total order of MPI collective sequences between processes, a check collective function CC is inserted before each MPI collective operation and before return statements. CC is depicted in Algorithm 2 in [START_REF] Saillard | Combining Static and Dynamic Validation of MPI Collective Communications[END_REF]. It takes as input the communicator com c related to the collective call c and a color i c specific to the type of collective. As multiple threads may call CC before return statements, this function is wrapped into a single pragma. Each function of a program is instrumented by Algorithm 3. If an error is about to occur, the program is stopped and an error message is returned with error type information.

Interprocedural Analysis

Because PARCOACH relies on an intraprocedural analysis, it miss errors across function boundaries and therefore it may produce false positive as well as false negative results. To extend PARCOACH with an interprocedural mechanism, we extended the intraprocedural approach through the application Call Graph (CG): nodes represent functions and edges model possible calls.

The main idea is to compute and reuse the summaries of each CFG through a CG traversal in reverse invocation order. For this purpose, the intraprocedural analysis is modified to return the valid sequence of collective operations for each function (ValidSeq). Based on this sequence of collective operations that all MPI processes ⊲ G: CFG, S, S ipw , S cc : sets created at compile-time 3:

if S ∪ S ipw ∪ S cc = / 0 then 4: ** STEP 1: Control flow errors detection ** 5:

for n in nodes containing a call to collective c do 6:

Insert call to CC(com c , i c ) before the call to c 7:

end for 8:

Before return statements insert 9:

# pragma omp single 10:

CC(communicator, 0)} 11:

* return O 12: end function will encounter, Algorithm 4 presents the interprocedural analysis. It takes as input the CG of a program and all the CFGs and returns the set O of conditionals that can lead to a deadlock.

Figure 5(a) shows the main CFG calling function g (whose CFG is depicted in Figure 5(b)). Performing our intraprocedural analysis would lead to a deadlock warning on main (for the MPI Barrier operations) and on g (for the MPI Allreduce collective). The resulting CFG after applying Algorithm 4 is illustrated in Figure 5(c): the call to g is replaced by the sequence of collective executed in g. Then, simply invoking the intraprocedural analysis on this new CFG results in no warning because there is one call to MPI Barrier on each path of the program. Therefore, the only warning issued by our combined analysis is related to the call to MPI Allreduce in function g (related to the if statement in node 2).

Experimental Results

We extended the PARCOACH implementation (GCC 4.7 plugin) to add analysis of hybrid applications. Thus, it is associated to GCC but simple to deploy in existing environments as it does not modify the compilation chain. To show the impact of PARCOACH analysis on the compilation and execution time, we tested the NAS-MZ [12], AMG benchmark [START_REF]CORAL Benchmarks[END_REF], the EPCC suite [START_REF] Bull | Performance Evaluation of Mixed-Mode OpenMP/MPI Implementations[END_REF] and HERA [START_REF] Jourdren | HERA: A hydrodynamic AMR Platform for Multi-Physics Simulations[END_REF]. All results were conducted on Tera100, a petaflopic supercomputer at the CEA. The interprocedural analysis has also been implemented and integrated in PAR-COACH (as a GCC plugin combined with a Python script). Figure 8 shows pieces of code from the IO module of the CORAL benchmark HACC. The calls in Figures 8(b) and 8(c) contain one and two calls to an MPI collective, respectively. Hence when the function in Figure 8(a) calls the others in different paths because of the switch, the execution could deadlock if the processes follow different paths. Of course the conditional statement does not depend on the rank number of each process, and therefore this is just a false positive. This interprocedural must be extended to a data flow analysis with the aim to study the dependence of these condition variables in order to know if they depend on the process rank or not.

Conclusion

The MPI+OpenMP approach is one solution to tackle the increasing node-level parallelism and the decreasing amount of memory per compute unit. Some production codes are already hybrid and other applications are in the development process. It is driven by available tools that could help debugging. That is why we developed the platform PARCOACH that helps application developers to check which interaction support is required for a specific hybrid code and checks the correct usage of blocking and non-blocking MPI collective communications in an MPI+OpenMP application. The main advantage of PARCOACH is that it highlights the statements responsible for the execution path potentially leading to future deadlocks or unspecified behaviors. We propose an adaptation of PARCOACH analyses to an interprocedural analysis. This enables us to reduce the number of false positives returned by the initial static analysis. However, this interprocedural analysis could be improved to propagate collective issue information and can be coupled to a data-flow analysis to avoid false positive results.

Fig. 3

 3 Fig. 3 Categories of possible errors in a hybrid program with N MPI processes and two threads per process.

Fig. 4

 4 Fig.4Automata of possible parallelism words. Nodes 0 and 2 correspond to code executed by the master thread or a single thread. Node 1 corresponds to code executed in a parallel region, and 3 to code executed in nested parallel region.

Algorithm 1

 1 Step 1: Static Pass of hybrid programs 1: function HYBRID STATIC PASS(G = (V, E), L) 2: ⊲ G: CFG, L: language of correct parallelism words 3: DFS(G, entry(G)) ⊲ parallelism words construction 4: MULTITHREADED REGIONS(G, L) ⊲ creates set S ipw 5: CONCURRENT CALLS(G) ⊲ creates set S cc 6: STATIC PASS(G) ⊲ creates set S 7: end function Static Pass Algorithm: To wrap-up all static algorithms, Algorithm 1 shows how analyses are combined. First the DFS function creates parallelism words.

Algorithm 3

 3 Step 2: Selective Static Instrumentation 1: function INSTRUMENTATION(communicator, G, S, S ipw , S cc ) 2:

Fig. 5

 5 Fig. 5 Example of interprocedural analysis

Fig. 6

 6 Fig. 6 Overhead of average compilation time with and without verification code generation

Fig. 7 Fig. 8

 78 Fig. 7 Execution-Time Overhead for MZ (NAS class B) and HERA with 8 threads per MPI process (Strong scaling)


Marwa Sridi
email: sridi.marwa@hotmail.fr

Bruno Raffin
email: bruno.raffin@inria.fr

Vincent Faucher
email: vincent.faucher@cea.fr

Cache Aware Dynamics Data Layout for Efficient Shared Memory Parallelisation of EUROPLEXUS

Keywords: EUROPLEXUS, Shared Memory, Cache-aware Data Layout, Parallel Programming 1

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

EUROPLEXUS (EPX)

is an industrial grade simulation software co-owned by the Commissariat à l' Énergie Atomique et aux Énergies Alternatives (CEA) and the Joint Research Center of the European Commission (EC/JRC). EPX si dedicated to the analysis of fast transient phenomena involving structures and fluids in interaction. It is based on a space discretization by means of Finite Elements, SPH Particles (Smooth Particle Hydrodynamics) or Discrete Elements for structures, or by means of Finite Elements, Finite Volumes or SPH Particles for fluids. Time integration is achieved through a conditionally stable explicit scheme. The solving algorithm is completely non-linear, at both geometric level (large displacements; large rotations) and material level (constitutive laws implementing plasticity or damage for example). The program provides a large number of kinematic links between entities, for instance for boundary conditions, contact between structures or fluid-structure interaction. EPX is characterized by its minimal use of non-physical parameters to enforce these links, such as penalty coefficients. It relies on direct methods to compute the link forces whenever it is possible and otherwise, the links are dualized by means of Lagrange Multipliers, the unknown forces being then deduced from the resolution of an additional linear system.

EPX is parallelized with MPI for clusters. But today supercomputer nodes show complex multi-core shared memory architectures. Relying on MPI for taking advantage of these multiple core lead to extra memory usage (phantom cells) and complex code to support load balancing. Task oriented programming models for shared memory architectures, available through programming environments like OpenMP, Intel TBB, Cilk or XKAAPI , offer an interesting alternative. The user only needs to delimit the potential parallelism through tasks or loop iterations, the runtime taking care of distributing these tasks to cores, relying on dynamics load balancing schedulers.

A large part of the computational cost in EPX classically occurs in the main loop iterating over the simulation elements to update their state (75% of the execution time for the MARA2 dataset introduced in Section 4). The cost of iterations is not uniform, mostly due to the different element formulations and materials co-existing within a single industrial model, but also to changes in the physical state of the dynamic systems modifying the work load required to compute the local equations. This calls for dynamics load balancing capable schedulers. Data access patterns also show some irregularities: updating the state of a given element often requires to access both its connected nodes and neighbor elements. The iterations of this main loop being independent, it can easily be parallelized with a task oriented environment.

The cost of accessing the data can significantly impair the performance. It is important to rely on data structures with a memory layout enabling to efficiently benefit from the cache hierarchy. However, the data layout EPX works on is inherited from upstream software tools used to mesh and partition the simulation domain. These tools can be different depending on the users and the type of simulation targeted. The quality of the data layout regarding the temporal and spatial locality is thus not ensured and can lead to numerous cache misses impairing performance.

In this paper we propose to have each task build its own local data structure gathering the data needed for processing the iteration range associated with this task. This lazy and parallel layout construction leads to very localized code modifications, and neither affects the global EPX data structure nor imposes constraints on the upstream tools for the simulation preparation. Experiments show that the cost of constructing these task local data structures is small compared to the compute cost and enables to significantly decrease the amount of cache misses. Globally our solution can lead up to a 40% performance gain when the task size is set to have the data fit in the L2 cache, compared to the simple loop parallelization on the original data structure. Without code modification we can benefit from the scheduler's dynamics load balancing capabilities. We show that our approach increases the benefits of load balancing by ensuring a strong locality of computation.

The remainder of the paper is organized as follows. In the second section we briefly present the EPX data structure. Then, we introduce our GROUP approach for both sequential and parallel implementation in section 3. In section 4 we provide some numerical experiments designed to evaluate the GROUP approach performance and we discuss the results. Section 5 concludes the paper and explores future research directions to further improve our implementation.

Dynamic Data Grouping

The EXP main loop iterates on the geometric elements that describe the simulation domain, to advance the simulation by one time step. Focusing on one element ei (Figure 1), carrying out an elementary computation from start to finish requires various data:

• ELEMENTARY data includes information allowing the definition of the mesh element ei through its intrinsic characteristics (deformation, internal variables. . . ).

• NODAL data includes information on the nodes of each mesh element (displacement, velocity, acceleration, forces. . . ).

• NEIGHBORHOOD data are variables that describe fluxes between the element and its neighbor cells (exchange between Ngi and ei in Figure 1), whether in form of energy or mass (for Eulerian or Arbitrary Lagrange Euler formulations). These concepts are valid in either structure only, fluid only or fluid/structure. These different data are stored in separate arrays. The order of the elements in these array is directly inherited from the files EPX loads for initializing the simulation. We have no guarantee that this layout will actually lead to an efficient use of the caching mechanism. The problem is amplified in a parallel context where false sharing for instance, i.e. the concurrent write/write or read/write of distinct data addresses but that belong to the same cache line, can lead to extra cache line invalidations.

To improve data locality and thus on cache efficiency, we propose to rewrite the main loop in 2 nested loops (Algo. 1). The outer loop iterates on GROUPS of elements (DO LOOP GR...) , while the inner loop (DO IELOOP...) iterates on each element of the processed loop. Each GROUP consists in a fixed number of elements.

This technique is thus based on a switch from the global structure of arrays layout (SOA) [START_REF] Govender | Development of a convex polyhedral discrete element simulation framework for {NVIDIA} kepler based {GPUs}[END_REF] to a local one. Therefore, an hybrid model, balancing the former layout at the scale of the different GROUPS and the array of structures layout at the scale of a single GROUP of elements, has been used. This hybrid SOAOS [START_REF] Wald | Fast construction of sah bvhs on the intel many integrated core (mic) architecture. Visualization and Computer Graphics[END_REF] model promotes the memory locality since the EPX code is based on a loop iterating over all the mesh elements.

To ensure that the inner loop can process the GROUP's elements with a high data locality, each outer loop iteration starts by filling a group data structure that consists in all elements belonging to these GROUPS (ELEMENTARY data) as well as all NODAL and NEIGHBOR-HOOD data that are required to update these elements (GROUP CONSTR routine). We consider that two elements are neighbors if they share at least one common side on the grid mesh. We determine the list of the cell neighbors by checking all the common sides with the adjacent cells. Then, we check if the global indexes belong to the list of the elements of the current GROUP. If this is not the case, it corresponds to extra data that does not belong to any other element of the GROUP and thus needs to be explicitly added to the GROUP. However, there is no need to copy the complete package for this added cell. We simply copy information ensuring its neighbor role (information coming from the NEIGHBORHOOD data), because it is already completely defined as a local element in another GROUP and therefore it is computed there.

Once the inner loop finished, the updated data are written back to the EPX global data structure (UPDATE routine). The elements in each GROUP correspond to consecutive indexes in the EPX global data structure. A local index is used to indicate the rank of each element relative to the other elements of the same GROUP. Beside reorganizing the main loop, we thus also need to change the arguments of the subroutines called in the inner loop to access to the GROUP's data instead of the global data structure. Elementary operations inside subroutines are kept intact.

1

Algorithm 1: The main iteration loop rewritten with 2 nested loops working on GROUPS

1 Subroutine COMPUTE(A) 2 ! Initialization step 3 NGR = N ELEM N B ELT GR 4 IF (N GR × N B ELT GR < NELEM) THEN 5 NGR = NGR+ 1 6 ENDIF 7 DO LOOP GR = 1, NGR 8 ISTART=1 + NB ELT GR 9 ×(LOOP GR -1)

Parallel Implementation

Our implementation relied on the XKAAPI task oriented programming environment 2 . XKAAPI is a C++ library providing a low overhead work stealing scheduler. Parallelizing a loop with independent iterations simply require to rely on the KAAPI FOREACH callback. The KAAPIF SET POLICY function enables to control the scheduling algorithm used for parallelizing the loop. The scheduler can be set to be static. In this case the user specifies an initial distribution pattern of the loop iterations through KAAPIF SET DISTIBUTION. A possibility is to split the iteration range in p sub ranges, where p is the number of cores involved. When starting the execution all cores are responsible to process one range. If the work stealing scheduler is used, the actual parallelism is extracted on demand according to the core availability, leading to a dynamics work load balancing. The scheduling algorithm works as follow: The first core initially gets the full iteration range and starts processing the loop sequentially. When a core becomes idle, it randomly selects a victim core. If this victim core is in charge of executing an iteration range, the idle core steals half of this iteration range. If the victim has no work to steal, an other victim is randomly selected. A grain parameter (KAAPIF SET GRAIN) prevents that an iteration range is split in two if too small, to avoid a performance degradation: the cost of stealing would not be compensated by the extracted parallelism.

The work stealing scheduling algorithm has proven performance [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF]. Today several parallel programming environments like Cilk [START_REF] Blumofe | Cilk: An efficient multithreaded runtime system[END_REF], TBB [START_REF] Reinders | Intel Threading Building Blocks[END_REF] or XKAAPI [START_REF] Gautier | KAAPI: A Thread Scheduling Runtime System for Data Flow Computations on Cluster of Multi-Processors[END_REF] rely on low overhead work stealing schedulers. The decentralized model of work stealing enable to ensure a good scalability at high core count in opposite to dynamics scheduler that rely on a centralized list cores refer to when they need more work. We made the choice to rely on XKAAPI, but the implementation could have been equally done with OpenMP. OpenMP offers all necessary constructs and several OpenMP runtimes rely on work stealing schedulers. In particular the XKAAPI scheduler can be used as runtime for OpenMP [START_REF] Broquedis | Libkomp, an efficient openmp runtime system for both fork-join and data flow paradigms[END_REF].

We first parallelized the original code (without group construction) with KAAPI FOREACH on the loop on the elements. For the version with group construction (Algo. 1), the loop on the groups (line 7) was parallelized with a KAAPI FOREACH (Algo. 2), the work inside a group being sequentially executed by the core in charge of that group (Algo. 2). Thus in parallel cores construct groups, compute the new states of the elements of their group accessing only data local to their group and write back the updated new element values to the global data structure. Because we expect that the work related to one group be sufficiently large, a grain of one can be used: a core can steal an iteration loop to a victim if this victim has at least 2 iterations left (the one currently processed and one waiting to be executed).

Experiments

Environment and Use Cases

All experiments are performed on a machine with 4 eight-cores CPU Intel Xeon clocked at 2.2 Ghz (hyperthreading disabled). Each core has its own L1 and L2 cache (256 KB) . The L3 cache (16MB) is shared among the eight cores of each CPU. We use XKAAPI version 3.0.6 and EUROPLEXUS development version. For all experiments, each thread is pinned on a different core of the machine. For this reason, in the following the term core is also used to designate the thread running on this core. All reported measures are the average over 5 runs to obtain significantly stable values. Experiments are based on the so called MARA2 simulation (Figure 2) named after an experiment in the field of nuclear safety and consisting in the relaxation of an explosive in a closed vessel. The fluid in the vessel is composed of a liquid covered with a gas, with a compressed bubble of another explosive gas located at the center of the vessel. The pressure

Results

We first evaluate the cost of building a group for different group sizes. Figure 3 plots the average time spent to build one group divided by the group size (Algo 1 routine GROUP CONSTR). These performance results are obtained for both, a sequential run (single thread) and one run parallelised on 32 threads. The results of both versions shows that the sequential version prevails slightly the 32 threads one. This is due to the intensification of the memory traffic when a large number of threads accesses simultaneously the shared memory. Notice that for the largest group size the performance for the sequential an parallel run are identical. In this case we only have one single group, and thus the parallel execution is reduced to a sequential one. This plot mainly highlights that the cost of building very small groups is high will it stabilises around 10ms/MB for larger groups.

We now compare the execution time of the Reference execution that corresponds to the original implementation using a single loop and no group, and the group execution that corresponds to the 2 nested loops with the group construction (Figure 4). In this last case the execution time includes the cost of the group construction as well as the update of the elements in the global EPX data structure. For group sizes above 0,25 MB, the GROUP version is about 35% faster than the reference one. The data rearrangement into each group leads to improved performance, mainly due to a more efficient cache usage. This also occurs for the largest group size where only one groupe is created. The variation of the GROUP version is closely related to the cost of group construction (Figure 3). Overall, the overheads of group creation and update of EPX main data structures stay small enough to keep significant performance gains.

We now run tests with 32 threads using both Mara mdm (Figure 5) and Mara big dataset (Figure 6). In each case we compare four different parallelizations: the Reference code (resp. GROUP code) with a static data partitioning (Reference-static, resp Group-static) and with work stealing enabled dynamic load balancing (Reference-steal, resp. Group-steal). The work stealing enabled Reference version outperforms the static one by about 20% for both, the Mara mdm and Mara big datasets. The GROUP code brings in most cases significant performance gains over the Reference code, the work stealing version often being the best one. Notice in particular that the performance of the static version can be affected by a load imbalance. For instance when there is only 63 groups (Figure 6) one thread has only one group while all the other have 2 groups each, causing a major performance drop. Work stealing enables to smooth this type of effect. It enables to take advantage of the uneven group processing time by allocating unprocessed groups to the fastest cores.

The acceleration with 32 threads is about 21 if we compare the performance for the Mara big dataset. It is 3 times greater than the reference acceleration. However, the speedup of the Mara mdm (Figure 4 and figure 5) simulation is 2 times lower than the Mara big one. This performance gap is due to the limited size of this dataset which reduces the parallelism that could be extracted in this case. We also looked at the average number of the L2 misses per core for the execution with 32 threads and groups of size 0.25 MB (L2 cache size) (Figure 7 and Figure 8). These L2 caches misses numbers were obtained using LIKWID3 [START_REF] Treibig | Likwid: A lightweight performance-oriented tool suite for x86 multicore environments[END_REF]. For each figure we provide the number for the execution of the full COMPUTE routine (Algo. 1) as well as the execution of the inner loop for the GROUP code (Algo. 1 DO IELOOP). The GROUP based approach enables to improve the number of L2 cache misses by about 15%, with about 35% of the remaining cache misses occurring in the inner loop, the others being spent when creating the groups (GROUP CONSTR routine) and updating the main data structure. (UPDATE routine). The work stealing based GROUP code generates slightly more cache misses than the static one. We attribute this to the coherency management over the L2 private cache and the L3 cache caused by the data migration during the steal operation.

Globally these results show that it pays off to dynamically build groups. The best results are obtained when the groups can tightly fit into the L2 cache. This lead to many small groups, enabling a fine grain load balancing but without the extra overheads that smaller groups lead to. The performance with larger group sizes is less stable and more dependent on the data set.

5 Related work M. Wolf [START_REF] Joseph | High Performance Compilers for Parallel Computing[END_REF] worked on algorithms to reorder nested loop to improve the cache efficiency by improving data locality. This approach requires changes in the code while working on the bounds of the loops as well as on their bodies. It is not always applicable especially in the case of codes with complex dependencies. M.Kandemir et al. [START_REF] Kandemir | Optimizing spatial locality in loop nests using linear algebra[END_REF] favored memory data reorganisation using linear algebra techniques to perform transformations on the matrices of data. P. Clauss [START_REF] Clauss | Parametric analysis of polyhedral iteration spaces[END_REF] proposed a more general method based on polynomial transformations applied to the data access. This requires the parameterisation of the loop bounds and the array sizes. The advantage of these approaches based on the reorganization of data tables is that the code modifications are restricted to the accessed data tables and do not affect the internal routines. However, they require an a priori knowledge on the order of data accesses.

Other researchers have worked on the data reordering in the memory to promote the temporal and the spatial locality. They apply specific algorithms of mesh elements renumbering like the Cuthill McKee [START_REF] Liu | Comparative analysis of the cuthillmckee and the reverse cuthillmckee ordering algorithms for sparse matrices[END_REF][START_REF] Cuthill | Reducing the bandwidth of sparse symmetric matrices[END_REF][START_REF] Cuthill | Several strategies for reducing the bandwidth of matrices[END_REF] algorithm used in some implementations of the mesh creator tool CAST3M. This algorithm permutes the rows of a sparse matrix to reduce the required bandwidth when accessing the data. The space filling curves (SFC) [START_REF] Sagan | Space-filling curves[END_REF] are another way to improve the cache use. They consist in projecting the multidimensional physical space in a onedimensional space, while maintaining as much as possible the spatial locality of elements. SFCs require to reorganize the global data structure with a priori knowledge on the access pattern. These methods are cache oblivious as the data layout is built independently of any specific cache parameter. Cache oblivious algorithms can be very efficient, but often outperformed by cache-aware approaches where the data layout is built knowing the cache parameters [START_REF] Yotov | An experimental comparison of cache-oblivious and cache-conscious programs[END_REF].

Conclusion

In this paper we propose to dynamically build groups gathering all the data needed for processing a given iteration sub-range. This approach enables to significantly improve performance while imposing only very local code modifications. We did not need to revisit the main data structures of the simulation software. Experiments show that our approach can lead to 40% performance increase on 32 cores when choosing groups that tightly fit in the L2 cache. We expect that this approach can be used in other simulation software than EPX on loops whose performance is impacted by excessive cache misses. Futur work will investigate the benefits of this approach for architectures with scratch pad. Next generation of Xeon Phi processor will have a NVRAM (the main memory) and a smaller high bandwidth RAM (the scratch pad). In this case the user needs to explicitly control the data to be stored in the RAM [START_REF] Michael | Two-level main memory co-design: Multi-threaded algorithmic primitives, analysis, and simulation[END_REF]. We expect that with minimal code modifications we can have our groups stored in the RAM while the globale data structure is stored in NVRAM, enabling to efficiently take advantage of the higher speed of the RAM.

Figure 1 :

 1 Figure 1: Representation of a mesh element ei and its neighbors in 2D

10 IEND= 11 CALL 15 CALL 5 NGR

 1011155 MIN(ISTART+NB ELT GR-1, NELEM) GROUP CONSTR (A → A ) 12 DO IELOOP = ISTART, IEND 13 CALL LOOPELM (A') 14 ENDDO UPDATE (A ← A) 16 ENDDO 17 END SUBROUTINE COMPUTE Algorithm 2: Skeleton of the XKAAPI code used to parallelize the outer loop on GROUPS (Algo. 1) 1 Subroutine COMPUTE(A) 2 ! Initialization step 3 NGR = N ELEM N B ELT GR 4 IF (N GR × N B ELT GR < NELEM) THEN

Figure 2 :

 2 Figure 2: Gas bubble expansion and pressure waves in the early moment of the MARA 2 simulation

Figure 3 :Figure 4 :

 34 Figure 3: GROUP construction costs for different sizes for sequential and 32 threads runs for Mara mdmFigure 4: Elapsed time for different group sizes for a sequential execution (Reference: single loop, no group; GROUP: 2 nested loops with GROUPS) for Mara mdm

Figure 5 :

 5 Figure 5: The elapsed time in the elementary loop with different GROUP sizes for a 32 threads run using the Mara mdm dataset

Figure 6 :

 6 Figure 6: The elapsed time in the elementary loop with different GROUP sizes for a 32 threads run using the Mara big dataset

Figure 7 :

 7 Figure 7: The L2 cache misses measurement within the EPX subroutines using the Mara mdm dataset

Figure 8 :

 8 Figure 8: The L2 cache misses measurement within the EPX subroutines using the Mara big dataset

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

https://code.google.com/p/likwid/


William L Griffin

Mathieu Chassé

N Nickel

Keywords: Ni Atomic number: 28 Atomic weight: 58.69 Isotopes radii: fourfold, 55 pm, fourfold square, 49 pm, fivefold, 0.63 pm, sixfold, 9 pm Pauling electronegativity: 1.91 First ionization energy: 737.1 kJ/mol Chondritic (CI) abundance: 1

come

Properties

Nickel has an atomic number of 28, an atomic weight of 58.69, and a density of 8.9 g/cm 3 . It is a silvery metal that melts at 1455 C and boils at 2913 C; it is magnetic, malleable, and ductile and conducts electricity (www.rsc.org/ periodic-table ).

Ni has five naturally occurring stable isotopes and up to 18 radioactive isotopes, all long extinct. Up to 0.4 permil difference in Ni-isotope composition is observed between coexisting kamacite and taenite in iron meteorites and is ascribed to diffusion during exsolution and cooling [START_REF] Lazar | Experimental determination of equilibrium nickel isotope fractionation between metal and silicate from 500 C to 950 C[END_REF]. However, the isotopic composition of Ni in chondritic meteorites and Earth is probably identical within current analytical capabilities [START_REF] Lazar | Experimental determination of equilibrium nickel isotope fractionation between metal and silicate from 500 C to 950 C[END_REF]. Significant fractionation of Ni isotopes can occur in low-temperature processes such as lateritic weathering of ultramafic rocks; progressive weathering leads to isotopically lighter soils (D 60 Ni Soil-Bedrock = À0.47 ‰) [START_REF] Ratié | Nickel isotope fractionation during tropical weathering of ultramafic rocks[END_REF].

Ni is a transition element, with the electronic configuration [Ar]3d 8 4s 2 ; consequently it shows both siderophile (ironloving) and lithophile (silicate-loving) behavior. In the crust and upper mantle, it normally occurs as Ni 2+ , whereas in the Earth's core, it is alloyed with iron.

History and Use

Nickel (Ni) was first isolated from NiAs (niccolite) and named by A. F. Cronstedt in 1751. The name derives from the Saxon "nickel," a malicious sprite ("Kupfernickel") who prevented them from extracting copper from niccolite, which they had mistaken for chalcopyrite (CuFeS 2 ). Its properties were determined by J.B. Richter in 1804.

The main use of Ni (ca 65 %) is for the production of stainless steel, with demand increasing by ca 5 %/annum (International Nickel Study Group; www.insg.org). It is also used in the production of nonferrous alloys or so-called "superalloys" which show better resistance to corrosion and are widely used in jet engines and turbines. Smaller proportions are used in various surface treatments such as electroplating and as hydrogenation catalysts.

Ores and Production

World usage of Ni amounts to between 1.8 and 2 million tonnes/year (2014)(2015).

Nickel is mined almost entirely from two types of deposit, magmatic sulfides and lateritic horizons, developed through the weathering of ultramafic rocks, usually in tropical areas (e.g., New Caledonia, Venezuela, the Dominican Republic, Brazil, Madagascar, and eastern Australia). Magmatic sulfides occur in ultramafic lava flows (komatiites; Figure 1) in Archean cratons, in large intrusions of layered gabbros, and in the Sudbury giant impact crater. In all of these, the major ore mineral is pentlandite (Ni,Fe) 9 S 8 , but Ni also forms a range of other sulfides, arsenides, and antimonides. While magmatic sulfide ores are important resources, >70 % of the world's supply of Ni is mined from lateritic deposits. During weathering, Ni concentrates in Ni-bearing limonite (Fe oxyhydroxide) and garnierite (referring to a group of green Ni-bearing magnesium phyllosilicates including serpentine, talc, sepiolite, smectite, and chlorite) in the middle- upper parts of the weathering profile, providing a resource that is easily mined and processed.

Geochemistry

The cosmic abundance of Ni estimated from carbonaceous chondrites is 10,910 ppm (Palme et al., 2014); estimates for its abundance in Earth lie between 17,000 and 19,000 ppm. Most of this is concentrated in the Earth's core; analyses of iron meteorites suggest the core contains ca 5 wt% Ni [START_REF] Mcdonough | 3.16 -Compositional model for the Earth's core A2[END_REF], leaving ca 1860 ppm in the mantle [START_REF] Palme | 3.1 -Cosmochemical estimates of mantle composition[END_REF] and 47 ppm in the continental crust [START_REF] Rudnick | 4.1 -Composition of the continental crust[END_REF]. Ni substitutes readily into the Mg-rich minerals of the mantle and is largely retained in the residue during partial melting. Ni concentrations in the crust are correspondingly quite low. The worldwide distribution of Ni in different rock types is shown in Figure 1; Figure 2 provides more detail on its distribution in mantle rocks. The variability shown in Figure 2 is easily understood in terms of the relative abundance of Ni in the important minerals of mantle rocks (Figure 3). This behavior is related to the charge and ionic radii of Ni. Under crustal and upper mantle conditions, the most common oxidation state of nickel is Ni(II). The effective ionic radius of Ni(II) depends on its coordination number ( [4] Ni, 0.55 Å; [4] Ni square, 0.49 Å; [5] Ni, 0.63 Å; and [6] Ni, 0.69 Å). The identical charge and similarity in ionic radii between [6] Ni and [6] Mg allow significant substitution of Mg by Ni in olivine, the main Ni bearer among mantle-related minerals (Figure 3). In addition, these distributions reflect mineral-mineral partitioning. Ni can substitute for Mg but not for Ca, which explains the higher Ni concentration in orthopyroxene compared to clinopyroxene (Figure 3). Ni substitution may also be favored in the Mg site of spinel, which is smaller than the Mg site of garnet.

The solubility of Ni in pure water is low (<2 mg/L), but in the presence of dissolved organic compounds, it can form complexes that are soluble. In anoxic conditions, it precipitates in sulfide phases (NiS, pyrite). Atmospheric concentrations of Ni are dominated by anthropogenic pollution, principally from the burning of oil, and have increased by a factor of >500 since the beginning of the Industrial Revolution [START_REF] Pacyna | An assessment of global and regional emissions of trace metals to the atmosphere from anthropogenic sources worldwide[END_REF]. Contributions to atmospheric Ni from smelting operations account for ca 10 % of the total; they have declined in recent years due to better emission controls, but previously devastated large areas in the vicinity of smelters at, e.g., Sudbury in Canada and the Kola Peninsula in arctic Russia.

Biological Utilization and Toxicity

Ni is known to be toxic in high doses, but they are reached only in special situations. However, Ni contamination of groundwater can represent a hazard, especially in acidic conditions where it can be leached from soils [START_REF] Nieminen | Biogeochemistry of nickel and its release into the environment[END_REF].

Nickel is an important micronutrient in the oceans, and dissolved Ni is variably enriched (ca 10 nM/kg) in deeper waters and depleted (ca 2 nM/kg) in surface waters, following the pattern of other micronutrients such as Fe and Zn, as well as P and Si. This distribution suggests that diatoms actively Nickel, Figure 3 Distribution of Ni concentration in various mantle-related minerals. (Data from GEOROC, Sarbas and Nohl (2008) and authors' unpublished data;105,630 analyses.) control the distribution of Ni in ocean waters [START_REF] Twining | Role of diatoms in nickel biogeochemistry in the ocean[END_REF]. The methanogenic bacteria that dominated Archean life appear to have been addicted to Ni. It has been proposed that a decline in the Ni contents of banded iron formations late in the Archean reflects a decrease in the availability of Ni in ocean water. This "nickel famine" may have led to the rise of the cyanobacteria, and their production of oxygen led in turn to the Great Oxidation Event at the end of the Archean [START_REF] Saito | Less nickel for more oxygen[END_REF][START_REF] Konhauser | Oceanic nickel depletion and a methanogen famine before the Great Oxidation Event[END_REF], allowing the evolution of multicellular life.

Summary

Nickel shows a very diverse range of geochemical behavior, being both a siderophile and a lithophile element. It has been important to mankind since the Bronze Age and is a mainstay of modern industrial society. It also is an essential element in the biogeochemistry of the oceans.

Cross-References

Figure 2

 2 Distribution of Ni concentration in various lithotypes (Data from GEOROC, Sarbas and Nohl (2008); 18,359 analyses). The boxes represent the first and third quartiles of the data distribution. The horizontal line represents the median, and the width of the adjacent notches shows the 95 % confidence interval about the median. The "whiskers" represent 1.5 x midrange (arithmetic mean of maximum and minimum values) beyond the first or third quartile Nickel, Distribution of Ni concentration in various mantle-related lithotypes. (Data from GEOROC, Sarbas and Nohl (2008); 5,403 analyses.)


Sofiane Bouzid

Magali Chaudey

Muriel Fadairo

Frédéric Perdreau

Strategic Groups in the French Franchising Sector

Keywords: Franchising, Strategic groups, Statistical models, Performance. JEL Classification: L21, C18

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


B Tanguy

J Besson

A Pineau

Comments on "Effect of carbide distribution on the fracture toughness in the transition temperature region of an SA 508 steel"

Keywords: A. Toughness, B. Statistical model, C. Finite elements analysis

 are recalculated using 2D and 3D FE simulations and new constitutive law, instead of small scale yielding assumption and likely improper choice for the behavior of SA 508 steel. A new model based on the weakest link concept with the determined carbide size distribution (CSD) is succesfully applied to predict K Jc measurements obtained on SA 508 steel and on a similar material (22NiMoCr3-7).

Introduction

This paper comments on a recent publication by Lee et al. [1] on cleavage fracture of SA 508 steel. Following Curry and Knott [2] who showed that the cleavage fracture toughness of a spheroidized steel depends on carbide distribution, most of the models for cleavage fracture of low alloy and mild steels are based on a statistical approach. Contrarily to mild steels, the full characterization of the carbide size distribution in quenched and tempered low alloy steels is difficult to obtain. This is why relevant experimental data are scarcely reported in the literature. In most of the carbide induced cleavage fracture models it is assumed that the carbide size distribution follows an "a priori" function (see e.g. [3,4]).

Recently in a comprehensive study by Lee and al. [1] (hereafter referred to as Lee) the carbide size distribution (CSD) of an SA 508 steel, which is a quenched and tempered low alloy steel, was given based on the analysis of more than 1400 particules. Lee has investigated the effect of carbide distribution on the fracture toughness, K Jc , of this material tested in the lower part of the transition temperature region where only cleavage fracture occurs without any prior ductile crack growth. Assuming that in this temperature range, cleavage fracture is controlled by a stress criterion and using the maximal principal stress profiles given by McMeeking [5] under plane strain small scale yielding (SSY) conditions, a linear relationship between the critical carbide size initiating cleavage fracture and the test temperature was found. Based on this and on the relationship between the critical nearest-neighbour distance and the CSD, a deterministic relation was proposed between K Jc and the CSD. However in Lee's work anomalously large values for the critical cleavage stress were reported (see Table 3 in [1]) considering the mechanical properties of this material compared to those obtained on the same steel [3] and in a recent study on a very similar steel [START_REF] Tanguy | ECF 14 , Fracture Mechanics beyond 2000[END_REF]. Moreover the SSY assumption for tests carried out on PCVN (Precracked Charpy V-Notch) specimens tested at fracture toughness values as large as 150MPa √ m (see Table 2 in [1]) appears to be strong.

Here, in order to comfort the results obtained by Lee [1], critical cleavage stress values have been recalculated based on the mechanical properties of SA 508 steel and using finite element (FE) calculations. Firstly the SSY assumption was kept in order to check the stress level values reported in Lee's work. Secondly, 3D numerical simulations of PCVN geometry were performed in order to validate the SSY assumption. Then, keeping the Lee's hypothesis that cleavage fracture in quenched and tempered bainitic steels is induced by carbides, and using the CSD obtained on SA 508 steel, a further extension of Lee's work is made to present a statistical model. This model follows the Beremin's formalism [3] but with the CSD given in Lee's work. Model parameters are adjusted to experimental results obtained on SA 508 steel with PCVN geometry and then applied to predict the K Jc values of another quenched and tempered bainitic steel, 22NiMoCr3-7 material which is very close to SA 508.

Materials and experiments

Table 1 gives the chemical compositions of the investigated materials. Both were quenched and tempered. Full details can be found in [1]. Details on experiments performed on SA 508 steel are given in [1]. Here it is simply reminded that elastic-plastic fracture toughness, K Jc , was determined using precracked Charpy V-notch (PCVN) specimens (10 × 10 × 55mm 3 ) (a 0 /W ∼ 0.5). Static ( ε = 10 -3 s -1 ) tensile tests were performed on smooth specimens in 22NiMoCr3-7 material to determine the stress-strain curves at various temperatures. Static plane strain fracture toughness tests were carried out on CT(1T) specimens according to standard ASTM E 1921-97 with a nominal crack length to specimen width ratio (a 0 /W ) of 0.5. These specimens tested between -90 • C and -30 • C led to pure cleavage fracture.

Results

Tensile properties

The evolution of the yield stress, σ Y , and ultimate tensile stress, R m , as a function of temperature is reported in fig. 1a for both steels. The results concerning SA 508 steel are taken from [1]. As shown in fig. 1, σ Y and R m are slightly lower for SA 508 steel than for 22NiMoCr3-7 steel. It is worth noting that at a given temperature, the difference between σ Y and R m is nearly the same between both steels, which indicates that the hardening capacity for both materials is quite similar over all the investigated temperature range. The strain hardening exponent, n, in ref. [1], was also given at different temperatures. However, the derivation of the function linking the flow stress and the plastic strain was not given in this paper. Assuming that the usual relation σ eq = K × ε n was used by the authors and using the Considère criteria, from which the relation R m = K(n/e) n is obtained, it is possible to determine the parameter K at each temperature (the same method was kept at -196 • C where fracture occurred probably before necking) using the tensile properties given for SA 508 steel. Corresponding stress-plastic strain curves obtained at -196 • C and -100 • C are plotted in fig. 1b (thick lines). These curves will be referred to as CL1 model in the following. On the same figure, the experimental stress-plastic strain curve for 22NiMoCr3-7 steel at -100 • C is also plotted. The stress-strain curves for both materials are very close. In the same figure, two other equivalent stress-plastic strain curves are also reported. These curves were obtained using the yield function given by McMeeking [5] (see eq. 1).

σ σ Y 1/N = σ σ Y + 3E 2(1 + ν) εp σ Y (1)
Here it is worth mentionning that the McMeeking's FE solutions were used by Lee to obtain the maximum principal stress, σ pI , at each test temperature, and that these values were used to calculate the critical carbide size. This later point will be discussed in the following. Keeping the n values given in ref. [1] (e.g. n = N in eq. 1), and assuming that the Kirchoff stress tensor can be approximated by the Cauchy stress tensor, leads to much higher stress levels (thin lines on fig. 1b referred to as CL2 model in the following) than stress-strain curves based on SA 508 tensile properties. In particular the ultimate tensile stress inferred from eq. 1 is strongly overestimated at both temperatures. For simplicity's sake, results at other experimental temperatures (-140 • C and -75 • C ) not shown here lead to the same conclusions were drawn. 3.2 Maximum principal stress, σ max pI , determination

R m SA 508 R m 22NiMoCr3-7 σ Y SA 508 σ Y 22NiMoCr3-7 R m σ Y Temperature, T ( • C ) σ Y , R m (MPa) 0 -

SA 508 steel [1]

In Lee's work the σ pI values at fracture (critical cleavage stress) was presumably obtained at each test temperature by using the McMeeking's FE solutions [5] with CL2 model. Based on the critical stress obtained for 22NiMoCr3-7 and A508 steels with CT(1T) geometry [START_REF] Tanguy | ECF 14 , Fracture Mechanics beyond 2000[END_REF], it appears that the stress values indicated in Lee's work are anomalously high, especially when considering that they were obtained with a PCVN geometry [START_REF] Joyce | [END_REF]. It is reminded that the results presented by McMeeking are based on a SSY assumption, plane strain analysis, and imposing an asymptotic dependence on mode I elastic crack-tip singular field (for more details see [5]). Similar calculations to those presented by McMeeking were made in our study but using the tensile stress-plastic strain curves inferred from the material data given by Lee (labelled CL1).

FE simulations were performed using software Zébulon [8]. Quadratic elements (8 nodes) with reduced integration were used. The boundary layer radii were modified in order to verify SSY conditions varying from 6 mm to 225 mm for K Jc values included between 36.3 and 150 MPa √ m. σ pI profiles corresponding to the experimental K Jc values given by [1] are reported in fig. 2a. This stress reaches a maximum at a distance, X c . In fig. 2b, the stress profiles corresponding to CL2 model are reported. The maximum values of σ pI for each experimental K Jc value are given in Table . 2 where the values reported by Lee are also given. From these results, two main conclusions can be drawn : (i) The maximum principal stress, σ max pI , values reported by Lee are much higher than those obtained from simulations with SSY assumption using the stress-strain curves of SA 508 steel (CL1 model). (ii) The σ max pI values reported by Lee are very close to those obtained from simulations with SSY assumption using CL2 model. Similar values are obtained, except at -196 • C were simulations give lower values. From these results it is inferred that the results presented by Lee were obtained using CL2 model as constitutive law for SA 508. Based on fig. 1b, it appears that this hypothesis leads to stress levels much larger than those expected from the material properties. It is then concluded that the σ max pI values reported in Lee's work were largely overestimated, and, that the critical carbide sizes calculated from these values were overly underestimated. Plane strain and full 3D simulations of PCVN specimens were also performed. σ max pI values are reported in Table . 2 while σ pI stress profiles obtained at each test temperatures are shown in fig. 3. It is observed that when the test temperature is higher than -100 • C , the SSY assumptions are no longer valid. This result is observed for both plane strain and 3D simulations. Therefore the σ max pI value is lower than that predicted from SSY assumption. Table 2 Comparison between maximum principal stress values, (σ max pI ), obtained with small scale yielding assumption (SSY) and full 3D simulations of PCVN tests using two different constitutive laws (CL1 and CL2).

-75 • C , K Ic = 150MPa √ m -100 • C , K Ic = 81MPa √ m -140 • C , K Ic = 71MPa √ m -196 • C , K Ic = 36.4MPa √ m x/(K Ic /σ Y ) 2
-75 • C , K Ic = 150MPa √ m -100 • C , K Ic = 81MPa √ m -140 • C , K Ic = 71MPa √ m -196 • C , K Ic = 36.4MPa √ m x/(K Ic /σ Y ) 2
σ max pI (MPa) -196 • C -140 • C -100 • C -75 • C SSY,
-75 • C , K Jc = 150MPa √ m -100 • C , K Jc = 81MPa √ m -140 • C , K Jc = 71MPa √ m -196 • C , K Jc = 36.4MPa √ m Plane strain x/(K Jc /σ Y ) 2 σ pI /σ Y 0.015 0.01 0.005 0 5.0 4.0 3.0 2.0 1.0 0.0 -75 • C , K Jc = 150MPa √ m -100 • C , K Jc = 81MPa √ m -140 • C , K Jc = 71MPa √ m -196 • C , K Jc = 36.4MPa √ m 3D x/(K Jc /σ Y ) 2

22NiMoCr3-7 steel

In order to determine the σ max pI values corresponding to the K Jc values obtained with this steel, 3D simulations of CT(1T) tests were performed. A good agreement between experimental and simulated load-CMOD curves were obtained at different temperatures. The calculated values of σ max pI are reported in Table 3.

Determination of the critical parameters

Based on the obtained values for σ max pI and keeping the Lee's assumptions according to which cleavage fracture is stress controlled by the propagation of a microcrack located through the carbides thickness, the critical carbide sizes was calculated by the authors [1] using a modified Griffith equation given by :

d c = 4Eγ p π(1 -ν 2 )(σ max pI ) 2 (2)
where γ p is the effective surface energy (equal to 7 J/m 2 in Lee's work), ν is the Poisson's ratio. The calculated results based on PCVN specimens simulations using the constitutive equation CL1 for SA 508 steel and experimental tensile stress-plastic strain for 22NiMoCr3-7 steel are reported in Table 3. σ max pI values obtained with both plane strain and 3D simulations are reported in Table . 3 for SA 508 steel. These results show that the σ max pI values obtained from 3D simulations are very close to those calculated from plane strain assumption. The evolution of the critical carbide sizes with temperature is reported in fig. 4 where Lee's results are also included. As already indicated by Lee for SA 508 steel, the critical carbide size, d c , increases linearly with increasing test temperature, but for SA 508 steel the values obtained in the present study are higher due to the lower σ max pI values obtained. Referring to the CSD given by Lee, at the highest test temperature, a very few number of the carbide particles will be eligible to participate to the nucleation of microcracks. Results obtained for 22NiMoCr3-7 steel tested at higher temperature show a weaker temperature dependence of the critical carbide size. However it is observed that at similar temperature (∼ -80 • C ), both materials lead to similar values for the critical size of carbides particles (∼ 0.7µm), as expected due to the similarity of these materials.

Statistical model

Lee showed that the probability of finding carbides larger than a given size, d, is given by :

T ( • C ) σ max pI (MPa) X c (mm) d c (
P (size > d) = exp - d -d u d 0 m ( 3
)
with d u =0.00917µm, the size of smallest carbides observed, d 0 =0.10158µm and m = 1.192, the shape factor of the Weibull distribution.

Assuming that a carbide can lead to the nucleation of a microcrack as soon as plasticity occurs [9], the microcrack size distribution can be represented by the CSD. Fracture occurs when, for a given local stress, σ pI , the carbides size is greater than a critical value, d c . The value of d c can be linked to a local fracture toughness, k Ic by the following formula :

d c ∼ 1 π k Ic σ pI 2 (4)
Using eq. 3, the corresponding failure probability is :

p r = exp - d c -d u d 0 m ( 5
)
Using a statistical analysis similar to that one proposed by Beremin [3], it can easily be shown that the failure probability is given by :

P R = 1 -exp - V V 0 exp - d c -d u d 0 m ( 6
)
(where V is the stressed volume) or expressed in terms of stress as :

P R = 1 -exp   - V V 0 exp   -   1 σ 2 pI -1 σ 2 u 1 σ 2 0   m     (7)
where

σ 0 ∼ 1 π 1/2 k Ic d 1/2 0
. In Eq. 6, the carbide size threshold, d u , is equivalent to the introduction of a stress threshold, σ u , such that

d u ∼ 1 π k Ic σu 2
. This expression leads to P R → 0 when the stress is close to 0 and to P R → e V/V 0 -1 e V/V 0 when σ pI → σ u . For sufficiently large values of V /V 0 , P R tends toward 1.

Adjustment of the model parameters

Among the four model parameters, σ u , σ 0 , m and V 0 , the m value is given by Lee, i.e. m = 1.192. It is assumed that this value represents also the CSD for the 22NiMoCr3-7 steel. The reference volume, V 0 is chosen following Beremin's work [3], i.e. as a cubic volume containing about 8 prior austenite grains : 50 × 50 × 50µm 3 . Parameters σ u and σ 0 were fitted to data obtained from PVCN tests on SA 508 steel. Actually only the σ u parameter has to be fitted on experimental data because it can easily be shown that σ u /σ 0 = (d 0 /d u ) 1/2 . In SA 508 steel, σ u /σ 0 = 3.328 was obtained from Lee. In order to investigate the predictive capability of the model the parameter σ u was determined using SA 508 fracture toughness measured at -100 • C such that a failure probability close to 50% (45%) was obtained for K Jc = 82.9MPa √ m.

Results and discussion

Figure 5a shows the evolution of the P R as a function of the

Summary -Conclusions

(1) Anomalously large values for the cleavage stress reported by Lee et al [1] are likely due to two effects in their analysis : (i) the use of a constitutive law for SA 508 steel in Mc Meeking numerical calculations which largely overestimates the stress-strain curves for this material; (ii) the use of plane strain SSY assumption which does not apply when the fracture toughness (i.e. test temperature) is too high. (2) Lee's results concerning carbide size distribution are used in a new statistical model based on the weakest link concept to predict the variation of fracture toughness with temperature in SA 508 steel.

(3) This model is also applied to another steel (22NIMoCr3-7) for which a larger data base was available, using the CSD measured by Lee. A good agreement between experimental and predicted scattering for the fracture toughness is obtained provided that the parameters appearing in this model are correctly fitted.

Fig. 1 .

 1 Fig. 1. Evolution of the tensile properties of SA508 and 22NiMoCr3-7 steels as a function of temperature. a) Yield stress and ultimate tensile stress b) Flow stress with CL1 and CL2 models.

Fig. 2 .

 2 Fig. 2. Maximal principal stress profiles obtained in SSY conditions for SA 508 steel a) with CL1 model b) with CL2 model

Fig. 3 .

 3 Fig. 3. Principal stress profiles in the mid-section of a PCVN specimen corresponding to the experimental fracture toughness values for SA 508 steel. a) Plane strain b) 3D simulations.

Fig. 4 .

 4 Fig. 4. Critical cardide size (d c ) vs test temperature for SA 508 steel and 22NiMoCr3-7 steel.

- 75 •KFig. 5 .

 755 Fig. 5. Evolution of the failure probabilities obtained with the statistical model. a) σ u = 14000 MPa fitted on SA 508 results at -100 • C on PCVN geometry (σ 0 = 4200MPa, n = 1.192, V 0 = 0.000125mm 3 ) b) Prediction of the fracture toughness scattering on 22NiMoCr3-7 steel with CT(1T) geometry. ρ is the initial mesh crack tip radius.

Table 3

 3 Maximal principal stress (σ max pI ), distance (X c ) from precrack tip to location of maximum principal stress, critical carbide size (d c ).

	µm)

K

  Jc obtained from the simulations of PCVN tests for SA 508 steel at the test temperatures. In this figure the experimental K Jc values are represented by a full circle. Using the K Jc value at -100 • C to fit σ u (fig.5a) leads to σ u = 14000MPa and σ 0 = 4300MPa. In fig.5a, it is shown that fitting the σ u parameter at -100 • C tends to high P R values for the experimental data obtained at lower temperatures. It is clear that testing the applicability of the present model to the results published by Lee[1] would require a set of experimental results much larger than that reported in their publication. Jc values corresponding to P R = 10%, 50% and 90% are shown in fig.5b. Using σ u = 14000MPa leads to a good prediction of the experimental scattering at all test temperatures investigated (fig.5b).

	The statistical model was then applied to predict the experimental toughness
	scattering obtained on 22NiMoCr3-7 steel for which a wider data base was
	available. The 3D numerical simulations of CT(1T) tests were post-processed
	in order to evaluate the failure probabilities. For each test temperature,
	the K