Ksander Ejjaaouani
email: ksander.ejjaaouani@inria.fr

Olivier Aumage
email: olivier.aumage@inria.fr

Julien Bigot
email: julien.bigot@cea.fr

Michel Méhrenberger

Hitoshi Murai
email: h-murai@riken.jp

Masahiro Nakao
email: masahiro.nakao@riken.jp

Mitsuhisa Sato
email: msato@riken.jp

Michel Mehrenberger

InKS, a Programming Model to Decouple Algorithm from Optimization in HPC Codes

Keywords: programming model, separation of concerns, HPC, DSL

Existing programming models tend to tightly interleave algorithm and optimization in HPC simulation codes. This requires scientists to become experts in both the simulated domain and the optimization process and makes the code difficult to maintain or port to new architectures. In this paper, we propose the InKS programming model that decouples these concerns with two distinct languages: InKS pia to express the simulation algorithm and InKS pso for optimizations. We define InKS pia and evaluate the feasibility of defining InKS pso with three test-languages: InKS o/C++ , InKS o/loop , InKS o/XMP . We evaluate the approach on synthetic benchmarks (NAS and heat equation) as well as on a more complex example (6D Vlasov-Poisson solver). Our evaluation demonstrates the soundness of the approach as it improves the separation of algorithmic and optimization concerns at no performance cost. We also identify a set of guidelines for the later full definition of the InKS pso language.

Introduction

It is more and more common to identify simulation as the third pillar of science together with theory and experimentation. Parallel computers provide the computing power required by the more demanding of these simulations. The complexity and heterogeneity of these architectures do however force scientists to write complex code (using vectorization, parallelization, accelerator specific languages, etc.). These optimizations heavily depend on the target machine, the code has to be adapted whenever it is ported to a new architecture.

As a result, scientists have to become experts in the art of computer optimizations in addition to their own domain of expertise. It is very difficult in practice to maintain a code targeting multiple distinct architectures. One fundamental cause for this situation is the tight interleaving of two distinct concerns imposed by most programming models. On the one hand, the algorithm comes from the expertise of the domain scientists and does not depend on the target architecture. On the other hand, optimizations form another domain of expertise and has to be adapted for a given architecture. Thereby, both algorithm and optimizations concerns are expressed within a single code. This mix impedes simulation codes maintainability and readability while hindering developers productivity.

Many approaches have been proposed to improve this situation in the form of libraries or languages [START_REF] Lee | Implementation and performance evaluation of Xcal-ableMP: A parallel programming language for distributed memory systems[END_REF]4,[START_REF] Tang | The Pochoir stencil compiler[END_REF][START_REF] Steuwer | LIFT: A functional dataparallel IR for high-performance gpu code generation[END_REF]. Approaches based on automated optimization processes typically isolate the algorithmic aspects well, but restrict their domain of applicability and/or the range of supported optimizations. Approaches based on optimization tools and libraries enable optimization specialists to express common optimizations efficiently but leave others mixed with the algorithm.

In this paper, we propose the Independent Kernel Scheduling (InKS) programming model to separate algorithm from optimization choices in HPC simulation codes. We define the InKS pia language used to express the algorithm of an application independently of its optimization. This separation aims to improve the readability and maintainability of codes while easing portability and new optimization expression. This approach is used for common optimizations while InKS o/C++ for less common optimizations. Such a program can then be optimized using InKS o/XMP and InKS o/loop , two domain-specific languages (DSLs) which ask for optimization information only. While these DSLs target some common optimizations, InKS o/C++ can be used for less common ones.

This paper makes the following contributions: 1) it defines the InKS programming model and its platform-independent algorithmic language InKS pia ; 2) it proposes an implementation of InKS and tests the InKS pso approach with three optimization DSLs, InKS o/C++ , InKS o/loop and InKS o/XMP ; and 3) it evaluates the approach on the synthetic NAS parallel benchmarks [START_REF] Bailey | The NAS parallel benchmarks[END_REF] and on the 6D Vlasov-Poisson solving with a semi-Lagrangian method.

The remaining of the paper is organized as follows. Section 2 analyzes related works. Section 3 describes InKS and its implementation. Section 4 shows the use of InKS on a 6D Vlasov-Poisson solver. Section 5 evaluates our approach. Section 6 concludes the paper.

Related Work

We now present approaches currently available to scientific programmers to help implementing simulation applications. A first widely used approach is based on imperative languages such as Fortran or C. Libraries like MPI extend this to distributed memory with message passing. Abstractions very close to the execution machine make fine-tuning possible to achieve good performance on any specific architecture. It does however require encoding complex optimizations directly in the code. As there is no language support to separate the algorithm and architecture-specific optimizations, tedious efforts have to be applied [START_REF] Höhnerbach | The vectorization of the Tersoff multi-body potential: An exercise in performance portability[END_REF] to support performance portability. Algorithm and optimizations are instead often tightly bound together in codes.

A second approach is offered by tools (libraries, frameworks or language extensions) that encode classical optimizations. OpenMP [START_REF] Chandra | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF], REPARA [START_REF] Danelutto | Introducing parallelism by using repara c++11 attributes[END_REF] or Kokkos [4] support common shared-memory parallelization patterns. For example, Kokkos offers multidimensional arrays and iterators for which efficient memory mappings and iteration orders are selected independently. UPC [START_REF] El-Ghazawi | UPC: Distributed Shared Memory Programming[END_REF] or XMP [START_REF] Lee | Implementation and performance evaluation of Xcal-ableMP: A parallel programming language for distributed memory systems[END_REF] support the partitioned global address space paradigm. For example, in XMP, directives describe array distribution and communications between nodes. These tools offer gains of productivity when the optimization patterns they offer fit the requirements. The separation of optimizations from the main code base also eases porting between architectures. Even if expressed more compactly optimizations do however remain mixed with the algorithm. For instance, in OpenMP or REPARA, parallel concerns are specified on top of an existing code which already carries optimization choices, such as loops order.

A third approach pushes this further with tools that automate the optimization process. For example, PaRSEC [START_REF] Hoque | Dynamic task discovery in PaRSEC: A data-flow task-based runtime[END_REF] or StarPU [START_REF] Augonnet | StarPU: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] support the manytasks paradigm. In StarPU, the user expresses its code as a directed acyclic graph (DAG) of tasks with data dependencies that is automatically scheduled at runtime depending on the available resources. Other examples are SkeTo [START_REF] Tanno | Parallel skeletons for variable-length lists in SkeTo skeleton library[END_REF] or Lift [START_REF] Steuwer | LIFT: A functional dataparallel IR for high-performance gpu code generation[END_REF] that offer algorithmic skeletons. Lift offers a limited set of parallel patterns whose combinations are automatically transformed by an optimizing compiler. Automating optimization improves productivity and clearly separate these optimizations which improves portability. The tools do however not cover the whole range of potential optimizations such as the choice of work granularity inside tasks in StarPU for example. The algorithm remains largely interleaved with optimization choices even with this approach.

A last approach is based on DSLs that restrict the developer to the expression of the algorithm only while optimizations are handled independently, such as Pochoir [START_REF] Tang | The Pochoir stencil compiler[END_REF] or PATUS [6], DSLs for stencil problems. In Pochoir, the user specifies a stencil (computation kernel and access pattern), boundary conditions and a space-time domain while all optimizations are handled by a compiler. These approaches ensure a very good separation of concerns. The narrower the target domain is, the more efficient domain and architecturespecific optimizations are possible. However, it makes it less likely for the tool to cover the needs of a whole application. On contrary, the wider the target domain is, the less efficient optimizations are possible.

To summarize, one can consider a continuum of approaches from very general approaches where the optimization process is manual to more and more domain specific where the optimization process can be automated. The more general approaches support a large range of optimizations and application domains but yield high implementation costs and low separation of concerns and portability. The more automated approaches reduce implementation costs and offer good separation of concerns and portability but restrain the range of supported domains and optimizations. Ideally, one would like to combine all these advantages: 1) the domain generality of imperative languages, 2) the ease of optimization offered by dedicated tools and 3) the separation of concerns and performance portability offered by DSLs with 4) the possibilities of fine and manual optimizations offered by both imperative languages and dedicated tools. The following section describes the InKS programming model that aims to combine these approaches to offer such a solution.

The InKS programming model

This section first introduces the design of our InKS programming model, based on the use of distinct languages to express the algorithm and optimization choices separately. It then present a prototype implementation of the model composed of InKS pia , the algorithm language, and two InKS pso optimization The simulation algorithm consists in the set of values computed, the formula used to produce them as well as the simulation inputs and outputs. Optimization choices include all that is not the algorithm: e.g. the computing unit selected for each computation, their ordering, the memory location for each value, etc. Multiple optimization choices can differ in performance but simulation results depend on the algorithm only. The InKS approach is summarized in Figure 1. The InKS pia language is used to express the algorithm with no concern for optimization choices.

A compiler can automatically generate nonoptimized choices from an InKS pia specification, mostly for testing purposes. The InKS pso language is used to define optimizations only, while other information is gathered from the InKS pia code.

Many versions of the optimization choices can be devised for a single algorithm; for example to optimize for multiple targets. We now describe the InKS pia language. and proposes three preliminary DSLs to validate strategies for the design of the InKS pso language: InKS o/XMP handles domain decomposition, InKS o/loop focuses on efficient loop nest and InKS o/C++ enables to use C++ and write arbitrary complex optimization.

3.1 The InKS pia language.

In InKS pia [START_REF] Aumage | InKS, a programming model to decouple performance from semantics in simulation codes[END_REF] values are stored in infinite multidimensional arrays based on dynamic single assignment (DSA, each coordinate can only be written once). Memory placement of each coordinate is left unspecified. Computations are specified by kernel procedures that 1) take as parameters data arrays and integer coordinates; 2) specify the coordinate they might read and will write in each array; and 3) define either a C++ or InKS implementation. An InKS implementation defines kernels validity domains: coordinates where C++ kernels can generate values in arrays. Kernel execution order is left unspecified. The simulation entry point is a kernel marked public. Listing 1 presents a simple InKS pia code: a 1D 3-point stencil computation. The simulation consists of one 2D logical array, Array (line 12), two kernels, stencil3 (line 1) and boundary (line 6), and is parameterized by two integers, X and T (line 11). Line 12 specifies that the simulation starts with a subset of Array (every values in the space dimension, at the first time-step) and expects, as output, another subset: every values in the space dimension at the last time-step. To do so, it can call the stencil3 and boundary kernels with a specific set of values for their integer parameters x and t (validity domain) and with Array as their logical array parameter, as expressed on line 15 and 16.

A InKS pia code specifies a parameterized task graph (PTG) [7]. This graph is encoded using the Polyhedron Model [START_REF] Feautrier | Polyhedron Model[END_REF]. It provides a compact representation of static control programs. This representation covers a large range of problems but imposes a few limitations. Mostly, all the problem parameters must be known at launch time, it does for example not support adaptive mesh or time-steps. It is still possible to express these concerns outside InKS pia and call the InKS implementation multiple times with different parameters.

Once the algorithm is specified in InKS pia , the goal is to write the optimization choices meaning choosing a memory layout and a scheduling of the kernels. Two approaches exist to produce these choices: the automatic compiler or InKS pso . All approaches rely on the Kernel file. This file is generated by the InKS pia compiler which translates the InKS pia kernels to C++ functions. The InKS pia compiler can produce generic optimization choices with a valid but non-optimized computations scheduling and memory allocations to execute them. Scheduling and memory layouts are computed using the Integer Set Library [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF] and recent works on modular mapping in the polyhedron model [START_REF] Isoard | Extending polyhedral techniques towards parallel specifications and approximations[END_REF]. Arbitrarily complex versions of optimization choices can also be written manually in plain C++. These functions can be called from any ex-

1 kernel stencil3(x, t) : ( 2 double A(2) {in: (x-1:x+1, t-1) | out: (x, t)} 3 ) 4 #CODE(C) A(x, t+1) = 0.5*A(x, t-1) + 0.25*(A(x-1, t-1)+A(x+1, t-1)); #END 5 6 kernel boundary(x, t) : ( 7 double A(2) {in: (x, t-1) | out (x, t)} 8 ) 9 #CODE(C) A(x,t) = A(x, t-1); #END public kernel inks_stencil(X, T) : ( double Array(2) {in: (0:X, 0) | out: (0:X, T-1)} ) #CODE(INKS) stencil3 (1:X-1, 1:T) : (Array), boundary {(0, 1:T), (X-1, 1:T)} : (Array) #END
Listing 1: 1D stencil computation on a 2D domain in InKS pia isting code whose language supports the C calling convention. However, that approach requires information present in InKS pia to be repeated. The InKS pso DSL thus interface the optimization process with InKS pia , offering the optimization specialists to specify optimization only. InKS o/XMP (illustrated in Listing 2) handles distributed memory domain decomposition by combining C and directives based on XMP and adapted for InKS. The compiler replaces these directives by C and XMP code. The inks decompose directive supports static or dynamic allocation of logical arrays described in the algorithm. The domain size is extracted from InKS pia source and the user only has to specify its mapping onto memory. As in XMP, InKS o/XMP supports domain decomposition mapped onto an XMP topology. In InKS pia code, there are no concerns for memory optimization such as dimension ordering or memory reuse. Therefore, InKS o/XMP supports dimension reordering and folding which consists in reusing the same memory address for subsequent indices in a given dimension. The exchange directive supports halo exchanges. The user specifies which dimension should be exchanged and which computational kernel will be executed after the exchange. From these informations and the InKS pia kernel specification, the InKS o/XMP compiler then computes the halo size. While XMP requires halo values to be stored contiguously with the domain, InKS o/XMP supports a dynamic halo extension where halo values are stored in dedicated, dynamically allocated buffers to reduce memory footprint.

InKS o/loop (illustrated in Listing 3) offers to specify manually loop nests for which the compiler generates plain C++ loops. Plain C++ is usable in combination with InKS o/loop . The loop keyword introduces a nest optimiza- Listing 3: InKS o/loop optimization choices of a loop from the InKS pia code presented in Figure 1 with its C++ result and use tion with a name, the list of parameters from the algorithm on which the loop bounds depend and a reference to the optimized kernel. Loop bounds can be automatically extracted from InKS pia , but the Set keyword makes it possible to restrict these bounds. The Order keyword specifies the iteration order on the dimensions named according to the InKS pia code. The Block keyword enables the user to implement blocking. It takes as parameters the size of block for the loops starting from the innermost one. If there are less block sizes than loops, the remaining loops are not blocked. The Buffer keyword supports copying data in a local buffer before computation and back after to ensure data continuity and improve vectorization. The compiler uses data dependencies from the InKS pia code to check the validity of the loops order and generate vectorization directives whenever possible.

             ∂f (t, x, v) ∂t + v.∇xf (t, x, v) -E(t, x).∇vf (t, x, v) = 0 -∆φ(t, x) = 1 -ρ(t, x) E(t, x) = -∇φ(t, x) ρ(t, x) = f (t, x, v)dv

The 6D Vlasov/Poisson problem

The 6D Vlasov-Poisson equation, presented in (1), describes the movement of particles in a plasma and the resulting electric field. We study its resolution for a single species on a 6D Cartesian mesh with periodic boundary conditions. We solve the Poisson part using a fast Fourier transform (FFT) and rely on a Strang splitting (order 2 in time) for the Vlasov part. This leads to 6 1D advections: 3 in space dimensions (x 1 , x 2 , x 3 ) and 3 in velocity dimensions (v 1 , v 2 , v 3 ). Each 1D advection relies on a Lagrange interpolation of degree 4. In the space dimensions, we use a semi-Lagrangian approach where the stencil is not applied around the destination point but at the foot of characteristics, only known at runtime, as described in more details in [18].

The main unknown is f (f6D in the code), the distribution function of particles in 6D phase space. Due to the Strang splitting, a first half time-step of advections is required after f6D initialization but before the main time-loop. These advections need the electric field E as input. E is obtained through the FFT-based Poisson solver that in turn needs the charge density ρ as input. ρ is computed by a reduction of f6D. The main time-loop is composed of 4 steps: advections in space dimensions, computation of the charge density (reduction) and electric field (Poisson solver) and advections in velocity dimensions. The algorithm of the simulation is presented in Figure 2.

The remaining of the section presents two optimizations implemented in the Selalib [START_REF] Kormann | Massively parallel semi-lagrangian solution of the 6d vlasov-poisson problem[END_REF] version of the 6D Vlasov-Poisson problem.

The 6D nature of f6D requires a lot of memory, but the regularity of the problem means it can be distributed in blocks with good load-balancing.

[m][l][k][j][i] with t6d[n][m][l][k][j][i] #pragma xmp template_fix t6d[N][M][L][K][J][I] f6d = (double(*)[M][L][K][J][I]) xmp_malloc(xmp_desc_of(f6d), N, M, L, K, J, I);
Listing 4: 6D decomp. & halo exchange in InKS o/XMP and resulting code Advections account for the main computational cost of the problem, up to for 95% of the sequential execution time. Six loops surround the stencil computation of each advection and in a naive version, the use of a modulo to handle periodicity and application along non-contiguous dimensions slow down the computation. To enable vectorization and improve cache use, we copy f6D elements into contiguous buffers along with the left and right halos. Advections are applied on these buffers before copying them back into f6D. Blocking further improves performance by copying multiple elements at a time. Listing 5 corresponds to the InKS o/loop implementation of these optimizations in a sequential version and presents the generated code. The Poisson solver 1 /*** advec.iks InKSpia algorithmic file ***/ 2 kernel advection3 (i, j, k, l, m, n, t, K, step) : ( 3 double f6d {in: (i, j, 0:K, l, m, n, t, step-1) |out: (i, j, k, l, m, n, t, step)}, 4 int disp {in: (n)}, 5 double coef{in: (n, 0:4)} 6 ) 7 #CODE (C) 8 f6d(i, j, k, l, m, n, t, step) = 9 coef(n, 0) * f6d(i, j, k-1+disp(n), l, m, n, t, step-1) + coef(n, 1) * f6d(i, j, k+0+disp(n), l, m, n, t, step-1) + coef(n, 2) * f6d(i, j, k+1+disp(n), l, m, n, t, step-1) + coef(n, 3) * f6d(i, j, k+2+disp(n), l, m, n, t, step-1); #END public kernel main_code(t, I, J, K, L, M, N, Niter) : ( double coef(2), int disp(1), double f6d(6) {in: (i, j, 0:K, l, m, n, t, step-1) |out: (i, j, k, l, m, n, t, step)} ) #CODE(INKS) advection3 (0:I, 0:J, 0:K, 0:L, 0:M, 0:N, 1 relies on a remapping scheme, where the domain decomposition is modified between each FFT execution.

Evaluation

This section evaluates the InKS model on the NAS benchmark, a simple stencil code and the 6D Vlasov-Poisson problem. We have implemented the algorithm of the following programs using InKS pia .

-4 sequential NAS kernels (IS, FT, EP and MG), C++ version [START_REF] Griebler | Efficient NAS benchmark kernels with c++ parallel programming[END_REF] as reference; finite difference 3D heat resolution (7-point stencil) ([15] as reference); -6D Vlasov-Poisson, using Fortran/MPI Selalib [START_REF] Kormann | Massively parallel semi-lagrangian solution of the 6d vlasov-poisson problem[END_REF] as reference.

For the 3d heat equation solver, two strategies were implemented: one uses double buffering (Heat/Buf) and the other implements a cache oblivious strategy (Heat/Obl). To evaluate the InKS programming model on Vasov-Poisson 6D, we have conduct three experiments. The first one is based on optimization choices written in plain C++ and cover the whole simulation. The NAS CG kernel relies on indirections not expressible in the polyhedron model of InKS pia . Its implementation would thus have to rely on a large C++ kernel whose optimization would be mixed with the algorithm. For the same reason, the NAS FT kernels was only partially implemented. InKS pia can however be used to express all other NAS kernels, the 3D heat equation solver as well as the 6D Vlasov-Poisson algorithm. Even if not as expressive as C or Fortran, InKS pia , through the polyhedron model, handles static controls programs. This covers the needs of a wide range of simulation domains and offers abstractions close to the execution machine rather than from a specific simulation domain. More specifically, it supports programs expressible as Parameterized Task Graphs [7]: a directed acyclic graph of tasks (here kernels) in a compact representation whose dependence and scheduling are parameterized by integers fixed at the InKS entry point execution. Among others, it can express computations such as FFTs or stencils with input coordinates unknown at compile-time, as in 6D Vlasov-Poisson.

InKS separates the specification of algorithm and optimization in distinct files. Multiple optimization strategies can be implemented for a single algorithm, as shown for the 3D heat equation where each relies on a specific memory layout and scheduling. Similarly, based on the InKS pia algorithm implemented for the IS, EP and MG NAS kernels, we have developed multiple optimization choices based on InKS o/C++ and either InKS o/loop or InKS o/XMP . For more complex cases such as the 6D Vlasov-Poisson problem, we were able to develop four optimization choices based on one InKS pia algorithm: (1) InKS o/XMP , (2) C with XMP, (3) InKS o/loop and (4) InKS o/C++ . This proves, to some extent, that the separation of concerns is respected.

Finding the right metric to evaluate the easiness of writing a code is a difficult question. As illustrated in Listings 1 and 5 however, algorithm expression in InKS pia is close to the most naive C implementation where loops are replaced by InKS validity domains with no worry for optimization. The polyhedron model used to represent the InKS pia code is compatible with a subset of the C language. Therefore, we could have choose the C as the algorithm language. However, the optimization process would have suffer from such a choice. Indeed, while InKS pia is designed to remove optimizations through the use of DSA and fine granularity kernels, C enables its users to reuse memory and specify a total scheduling. Although it is possible to analyze the code and to retrieve the DSA form and the partial order of the code using the polyhedron model, it will be complicated for optimization specialists to find which loops can be broke or which arrays can be expended from a C code.

Concerning the specification of optimization choices, it is close to their expression in C++. Table 2 compares the GNU complexity score of InKS optimizations to the reference code. InKS scores are slightly better, which indicates that our language is not more complex than C++. The difference comes from the extraction of the computational kernels, placed in the algorithm, which hides parts of the complexity. In addition, the use of InKS o/C++ to write optimizations let optimization specialists reuse their preexisting knowledge of this language. Similarly writing the InKS o/C++ version of optimization choices for the InKS pia version of Vlasov-Poisson 6D is close to the expression to the same optimization in Fortran, in the reference version. These considerations should not hide the fact that some information has to be specified both in the InKS pia and InKS o/C++ files with this approach leading to more code overall.

On contrary, InKS o/XMP and InKS o/loop enables the developer to specify optimization choices only while algorithmic information is extracted from InKS pia code. This is illustrated by Listing 2 presenting the InKS o/XMP 6D domain decomposition and the XMP result. Both are equivalent, but the InKS o/XMP expects only optimization choices parameters. Hence, one can test another memory layout, such as a different dimension ordering, by changing only a few parameters, while multiple directives must be modified in XMP. Moreover, using the domain decomposition directive offered by InKS o/XMP and XMP code, one can derive a simple parallel code from a InKS pia algorithm, as we did with the 3D Heat solver. Note that since InKS o/XMP is a wrapper for XMP using InKS pia to retrieve some information, it is usable in any codes that can be expressed using InKS pia and optimized with XMP. Similarly, with InKS o/loop (Listing 5), developers can easily test different optimization choices that would be tedious in plain C++. This is what we did with the 3D Heat solver based on the double buffering technique. As shown on Figure 1, using InKS o/loop , we have implemented 3 versions of the loops: the same as reference, one based on a 2D cache-blocking and a last one using a 3D cache-blocking by modifying almost nothing in the InKS o/loop code. Since InKS o/XMP and InKS o/loop are respectively usable with C and C++, InKS does not restrict the expressible optimization choices: one can still implement optimizations not handled by our DSLs in C/C++. Moreover, operations such as halo size computation or vectorization capabilities detection are automatized using the InKS pia code. In summary, the approach enables optimization specialists to focus on their specialty which make the development easier.

Regarding performance, the InKS approach makes it possible to express optimizations that do not change the algorithm. Optimizations of the four NAS parallel benchmarks and 3D heat equation solver in InKS were trivial to implement and their performance match or improve upon the reference as presented in Table 2. Investigation have shown that Intel ICC 18 does not vectorize properly the reference versions of MG and FT. The use of the Intel ivdep directive as done on the InKS versions leads to slightly better performance.

For the full Vlasov-Poisson 6D problem, InKS pia and C++ enable us to implement the same complex optimization strategies written in the reference version. Therefore, our implementation match the reference in terms of performance, as shown on Table 2. For the v 1 advection, both the InKS o/C++ and InKS o/loop optimizations of the InKS code achieve performance similar to the reference as shown in Table 2. For the parallel aspects, the InKS o/XMP optimization offers performance similar to XMP as shown on Figure 4. MPI is faster on all cases compared to both XMP and InKS o/XMP . At the moment, it seems that XMP does not optimize local copies which slows down the Poisson solver. Besides, some XMP directives are based on MPI RMA which make the comparison with MPI Send/Receive complex. Still, MPI is much harder to program: more than 350 lines of MPI and Fortran are required to handle domain decomposition, remapping for FFT and halo exchange in Selalib while 50 lines in XMP and 15 in InKS o/XMP . However, InKS pia , InKS o/loop and InKS o/XMP have limitations. As mentioned earlier, InKS pia can not be used to express non static control program, such as the CG or the full FT NAS kernel. It also makes more complex the optimization of arrays passed in parameters (as input or output). InKS o/loop offers very limited option, making it unusable with complex loops such as the one in the 3D heat solver optimized using the cache-oblivious strategy. Similarly, loop nests calling multiple InKS pia computational kernels can not be expressed Although we want to address some issues in the InKS pia language, InKS o/loop and InKS o/XMP were preliminary tests before a real InKS pso implementation. The goal was to propose a way to express optimizations and retrieving the algorithmic information in InKS pia code. These two DSLs enables us to highlight a set of guidelines for the full definition of the InKS pso language. First, this definition must be based on the concepts that InKS pso must express, i.e. the optimizations related to memory (allocations and layouts) and to computations scheduling. With InKS o/XMP we tried some tests on the memory aspects, providing a directive to allocate a logical array and to reorder its dimensions. InKS o/loop focus on computations scheduling by adding strong constraint on the order of the computations. Secondly, in order to express these optimization concepts only, it must be bound to its algorithmic counterpart, express in InKS pia . In InKS o/XMP , the domain decompistion directive makes reference directly to logical arrays described in InKS pia , making possible the computation of the size of each dimension and the halo size depending of the data being accessed. In InKS o/loop , the Order keyword refers to the structuring variable of a computational kernel defined in InKS pia . Finally, we must work on the mean to express these optimization concepts with the most generality. For this part, we think that working with existing tools is probably mandatory for most complex strategy, such as XMP for the domain decomposition. Using another existing tool for computations scheduling, making it more general. Similarly, a declarative language such as InKS o/loop may not be the best strategy: to handle the interaction between memory and computation, an imperative language could be easier.

Conclusion and future works

In this paper, we have presented the InKS programming model to separate the algorithm and the optimization choices and its implementation supporting two DSLs: InKS o/loop for loop optimizations and InKS o/XMP for domain decomposition. We have evaluated InKS on synthetic benchmarks and on the Vlasov-Poisson solving. We have demonstrated its generality and advantages in terms of separation of concerns to improve maintainability and portability while offering performance on par with existing approaches.

While this paper demonstrates the interest of the InKS model, it still requires some work to further develop it. We will improve the optimization DSLs; base InKS o/loop on existing tools and ensure interactions with InKS o/XMP . This will be done within the scope of the development of a complete InKS pso DSL enabling its users to manage memory placement and computations scheduling. This planed DSL will truly separates algorithm from optimization choices in static programs while its compiler could offer static analysis of the program to ensure correctness. We also want to target different architectures to demonstrate the portability gains of the InKS model.

Fig. 1 :

 1 Fig.1:The InKS model languages. The simulation algorithm consists in the set of values computed, the formula used to produce them as well as the simulation inputs and outputs. Optimization choices include all that is not the algorithm: e.g. the computing unit selected for each computation, their ordering, the memory location for each value, etc. Multiple optimization choices can differ in performance but simulation results depend on the algorithm only. The InKS approach is summarized in Figure1. The InKS pia language is used to express the algorithm with no concern for optimization choices. A compiler can automatically generate nonoptimized choices from an InKS pia specification, mostly for testing purposes. The InKS pso language is used to define optimizations only, while other information is gathered from the InKS pia code.

3. 2

 2 InKS pso DSLs implementation Implementing a complete optimization DSL for the InKS model is a long-term objective. In this study, we present the two DSLs InKS o/XMP and InKS o/loop . Both DSLs solely describe optimizations, and get missing information from InKS pia code.

Fig. 2 :

 2 Fig. 2: The 6D Vlasov-Poisson algorithm

Fig. 3 :

 3 Fig. 3: Dynamic halo exchange representation on a 2D domain Halos are required to hold values of neighbors for the advections. Connected halo zones would increase the number of points in all dimensions and consume too much memory. Split advections mean that halos are required in a single dimension at a time though. We therefore use dynamic halos composed of two buffers, one for each boundary of the advected dimension (denoted "right" and "left").Figure 3 shows this optimization on a 2D domain. Listing 4 shows the InKS o/XMP implementation of this strategy on 6D Vlasov-Poisson.

Figure 3

 3 shows this optimization on a 2D domain. Listing 4 shows the InKS o/XMP implementation of this strategy on 6D Vlasov-Poisson.

  InKS o/XMP optimization choices of the InKS pia code presented in Figure1and the C + XMP result

	1 /***InKSo/XMP specification***/
	2 #pragma xmp nodes p[nNodes]
	3 void inks_stencil(T& Input, size_t X, size_t T){
	4	#pragma inks decompose % Array // allocation of "Array" algorithmic array
	5	(2*2, 1) // dimension reordering, time dimension is fold by 2 and not distributed
	6	with t onto p // block decomposition mapped on the XMP topology
	7
	8	for(int t=1; t<T; t++){
	9	boundary(Array, 0, t);
		for(int x=1; x<X-1; x++)
		stencil3(Array, x, t);
		boundary(Array, X-1, t);
		}
	}
	/***XMP + C result***/
	#pragma xmp nodes p[nNodes]
	#pragma xmp template t[:][:] // xmp 1d logical array
	#pragma xmp distribute t[*][block] onto p
	void inks_stencil(T& Input, size_t X, size_t T){
		#pragma xmp align Array[t][x] with t[t][x]
		#pragma xmp template_fix t[2][X]
		Array = (double(*)[X]) xmp_malloc(xmp_desc_of(f6d), 2, X);
		for(int t=1; t<T; t++){
		boundary(Array, 0, t);
		for(int x=1; x<X-1; x++)
		stencil3(Array, x, t);
		boundary(Array, X-1, t);
		}
	}
	Listing 2: 1 /*** stencil.iloop InKSo/Loop optimization choices file ***/
	2 loop stencil3_loop(t, X) : stencil3 { // set "t" value
	3	// "Set" not specified -> loop bounds are computed, with a fixed "t"
	4	Order: x; // order of the loop
	5 }
	6 /*** C++ result ***/
	7 void stencil3_loop(T& Array, size_t t, size_t X){
	8	for (int x=0; x<X; x++)
	9	stencil3(Array, x, t);
	}
	/*** InKSo/Loop result used ***/
	void inks_stencil(T& Input, size_t X, size_t T){
		for(int t=1; t<T; t++){
		boundary(Array, 0, t);
		stencil3_loop(Array, t, X);
		boundary(Array, X-1, t);
		}
	}

  Listing 5: A loop nest in InKS pia optimized in InKS o/loop

			:Niter, K, 2) : (f6d, disp, coef)
	/* ... */
	#END
	/*** advec.iloop InKSo/Loop optimization choices file ***/
	loop advection3_loops(t, I, J, K, L, M, N, Niter) : advection3 { // set "t" value
	// "Set" not specified -> loop bounds are computed, with a fixed "t"
	Order: n, m, l, k, i, j; // order of the loop
	Block: 16;	// blocking on the inner dimension j
	Buffer: f6d(3);	// copy the third dimension of f6d to a 1d buffer
	}

Table 1 :

 1 Then, we evaluate InKS o/XMP and InKS o/loop on Vlasov-Poisson separately as they target different optimizations and are not usable together currently. A first experiment focuses on the sequential aspects with the intra-node optimization of the v 1 advection using either InKS o/loop or InKS o/C++ . A second experiment focuses on the parallel aspects with the charge density computation, the Poisson solver and a halo exchange optimized either with C/XMP or with InKS o/XMP . Summary of the InKS o/C++ , InKS o/loop and InKS o/XMP optimization choices implemented based on six InKS pia algorithm: the 3D heat solver, 4 NAS kernels (EP, FT, MG, IS) and the Vlasov-Poisson 6D solver.

		Algorithm	InKS o/C++		InKS o/loop	InKS o/XMP
	Heat	Buf Obl	3 versions: Same as reference, 2D and 3D cache blocking
		EP
	NAS	FT MG
		IS
		Initialization
	VP6D	Advections Reduction		adv_v1
		Poisson solver
			Execution time (second)	Complexity
	Benchmark	Reference	InKS o/C++	Rel. dev.	Ref.	InKS
	NAS/FT	49.16 (±0.16)	44.00 (±0.05)	11.66%	6	5
		NAS/IS	1.91 (±0.00)	1.90 (±0.01)	0.53%	55	52
	NAS/MG	4.61 (±0.02)	4.29 (±0.01)	7.46%	20	12
	NAS/EP	53.52 (±0.02)	55.35 (±0.37)	-3.30%	19	19
	Heat/Buf	2.76 (±2.78)	2.93 (±3.53)	-6.22%	5	3
	Heat/Obl	1.20 (±1.16)	1.18 (±0.96)	1.81%	22	13
		VP6D	25.92 (±0.21)	21.92 (±0.16)	18.28%	N/A	N/A
	Benchmark	Reference (C++)	InKS o/loop	Rel. dev.	Reference (Fortran)
		NAS/IS	1.91 (±0.00)	1.91 (±0.01)	0.00%	N/A
	NAS/MG	4.61 (±0.02)	4.32 (±0.02)	6.84%	N/A
		VP6D	4.18 (±0.05)	4.23 (±0.06)	-1.10%	3.63 (±0.09)
	Benchmark	Reference	InKS o/XMP	Efficiency (# cores)
	NAS/EP	53.52 (±0.02)	14.00 (±0.13)	70.5% (4 cores)

Table 2 :

 2 Execution time of the InKS o/C++ , InKS o/loop and InKS o/XMP implementations of the sequential NAS benchmark, class B -Time/iteration of the 3D heat equation (7point stencil), size (1024 3 ) -Time/iteration of the InKS o/C++ and Fortran implementation of the sequential 6D Vlasov-Poisson, size (32 6 ). Median and standard deviation of 10 executions -GNU Complexity score of the implementation.

The table 1 summarizes the optimization choices we have implemented. All codes are compiled with Intel 18 compiler (-O3 -xHost), Intel MPI 2018 and executed on the Irene cluster (TGCC, France) equipped with 192 GB RAM, two Skylake 8168 CPUs per node and a EDR InfiniBand interconnect. The table

2

summarizes the optimization choices we have implemented.

  Weak and strong scaling for 3 parts of the Vlasov-Poisson solver up to 64 nodes (1 process/node) on a 32 6 grid divided among processes (strong scaling) or 16 6 grid per process (weak scaling). Median of 10 executions.using InKS o/loop . Moreover, it can not respects the constraints imposed by a specific memory mapping and offers only a few optimization. Adding new optimization strategies would require to add new keywords, such as unrolling, loop fusion, etc. InKS o/XMP offers memory allocations but no controls on the access to this memory, letting optimization specialists in charge of the good use of this memory. Besides, it is usable only with XMP.

	1.3										1
	1.2
	1
							1
	0.8
	0.6
							0.316				0.1
			XMP weak		InKSo/XMP weak
			MPI weak		XMP strong
		InKSo/XMP strong		MPI strong							Local copies		Remote copies
	1	2	4	8	16	32	64	8	16	32	64	1	2	4	8	16	32	64
		(a) f6D reduction				(b) Poisson solver	(c) Dynamic halo exchanges
	Fig. 4:


Djamel Eddin Abdelli
email: abdelli22djamel@yahoo.com

Thanh Hung Nguyen

Stéphane Clenet

Ahmed Cheriet

T T Nguyen

S Clénet

Stochastic Metamodel for Probability of Detection Estimation of Eddy-Current Testing Problem in Random Geometric

Keywords: Eddy Current Testing, Finite Volume Method, Stochastic Model, Least Angle Regression Method, Probability Of Detection

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

owadays, non-destructive testing (NDT) is an essential element of component quality qualification. For inspecting materials and components, several methods are developed. One of these methods is the Eddy Current Method (ECT). The qualification of this process is then required. The size of defects that can be detected should be determined. In practice, it appears that imperfections on the sensors or uncertainties on the material characteristics to be inspected modify the response of the NDT device. The measurement is not deterministic and varies around a targeted value. For this reason, a statistical study must be carried out to calculate the possibility of detecting (or not) defects under different operation conditions [START_REF] Deng | Model based POD techniques for enhancing reliability of steam generator tube inspection[END_REF][START_REF] Gratiet | Model assisted probability of detection curves: New statistical tools and progressive methodology[END_REF].

A Probability Of Detection (POD) can be estimated to quantify this capability of the NDT device to detect a defect.

The ECT system can be represented by electromagnetic equations such as Maxwell's equations. The difficulty of solving these equations analytically leads to use a numerical method to construct an accurate model. In this work, the Finite Volume (FV) method has been applied. As long as the input parameters are no longer nominal values, so the numerical model FV can be represented by a stochastic system [START_REF] Clénet | Uncertainty Quantification In Computational Electromagnetics[END_REF] with input parameter which are random.

The response of ECT system is no longer deterministic and a sampling technic like the Monte Carlo Simulation Method can be used to characterize the output and particularly to calculate the POD. A high number of calls of the numerical model FV are required. So repeating the model FV for a fairly large number of realizations is time consuming, especially when the geometry is modified because the mesh should be modified. To avoid the remeshing (RM), a geometric transformation (GT) method [START_REF] Mac | Résolution numérique en électromagnétisme statique de problèmes aux incertitudes géométriques par la méthode de transformation: Application aux machines électriques[END_REF] has been used which consists in changing the coordinates of the nodes without changing their connectivity. The GT method has been used to solve an electrokinetic problem [START_REF] Mac | Solution of Static Field Problems With Random Domains[END_REF] and more recently a magnetoelectric problem [START_REF] Nguyen | Influence of Material and Geometric Parameters on the Sensor Based on Active Materials[END_REF]. To approximate the response of the stochastic system, a stochastic metamodel [START_REF] Gratiet | Model assisted probability of detection curves: New statistical tools and progressive methodology[END_REF][START_REF] Clénet | Uncertainty Quantification In Computational Electromagnetics[END_REF][START_REF] Nguyen | Uncertainty Quantification Using Sparse Approximation for Models With a High Number of Parameters, Application to a Magnetoelectric Sensor[END_REF] is constructed based on a polynomial chaos expansion and the Least Angle Regression (LAR) method.

In this work the numerical POD of the defect depth is estimated using the Hit/Miss method in both methods GT and RM. To reduce the time calculation a stochastic metamodel is used in both methods (GT and RM) to approximate the responses. The two approaches are compared on a NDT stochastic example.

II. ECT AND FV MODEL

The ECT problem is modeled by the axisymmetric formulation of the magnetic vector potential " A " in the domain Ω is used as follows:

1 .

, in

j w σ µ   -∇ ∇ + = Ω     s A A J (1)
with w , s J , σ , and µ are the angular frequency, source current density, electrical conductivity and magnetic permeability respectively.

N tt rci wc tc σ µ a Fig. 1 Coil-tube geometry [START_REF] Song | Eddy current flaw characterization in tubes by neural networks and finite element modelling[END_REF], with rci, wc, tc are the inner radius, the width, height of the coil, tt the thickness of the tube, and the defect a (shaded area).

This ECT problem consists of a coil located in a steam generator tube with an external circumferential defect. The geometry is shown in Fig. 1.

To solve Eq.1 with the FV model, we subdivide the domain into a large number of elements (triangle or quadrilateral). The integration of Eq.1 over each element (P) of the domain Ω is required as it is shows in Fig. 2.

The integration of first term of Eq.1 over the element P can be done by using the divergence theorem as in:

1 1 .( ). ( ). . 1 , i i i i s N P i i i f d c η η µ µ η µ ∆ ∇ ∇ = ∇ - = ∆ ∆ ∫∫ ∫ ∑ r  A ds A n A A (2)
where i N A , P A are the nodal magnetic vector potential of the nodes i N (i=1:4) , P respectively, i f are the faces with r n its outward normal vector , i c ∆ are the distances between the centres of i N and P , i η ∆ are the lengths of the edges of the element. The non-orthogonal term is considered and interpolated as in [START_REF] Versteeg | An Introduction to Computational Fluid Dynamic[END_REF]. We first solve the Eq.1 for the unknown potential A for all elements of the mesh and then we calculate the impedance of coil which is the quantity of interest to detect the default. It can be calculated using Faraday's law and Stokes theorem [START_REF] Ida | electromagnetic and calculation of field[END_REF] as follow:

. .

j w Z d Ω = - Ω ∫ ∫  s A dl J (3)

III. PROBABILITY OF DETECTION ESTIMATION

A. Empirical POD

The determination of the empirical POD requires considerable experimental measurements. A large number of samples are required to have a representative statistical population for signal of the defect. This method is time consuming and very expensive [START_REF]Update on the Tools for Integrity Assessment Project[END_REF]. In this paper, we work with the numerical POD.

B. Numerical POD • Monte Carlo Simulation Method (MCSM)

The Monte Carlo Simulation method [START_REF] Clénet | Uncertainty Quantification In Computational Electromagnetics[END_REF][START_REF] Ghanem | Stochastic Finite Elements. A Spectral Approach[END_REF] consists in generating a sample of size which corresponds to M realizations of the input parameters, to solve M times the model and then in a postprocessing to give an estimation of the POD.

Its accuracy depends on sample size, and convergence is relatively low. Indeed, thousands of random realizations are often required to obtain a desired accuracy of the POD, which necessarily implies a high calculation time.

For the implementation of the MCSM, it is assumed that the input parameters are independent and have a uniform distribution.

• Stochastic model

As long as the calculated response of the ECT system is no longer nominal (deterministic), the notation of a random variable ( Z ) is used. So the stochastic model is given by the function:

( ) , a . Z f X = ( 4
)
where a is the depth of the defect (in the case without defect a 0 = and Z is denoted 0 Z ) and X the vector a finite number of random input parameters. To alleviate the notations, we replace in the following ( ) . The depths of the defect are 0.05, 0.4970 and 0.9675 mm of the thickness tube with a width of 1.5 mm. The defect is circular and coaxial with the coil.

, a Z f X = by ( ) Z X (or Z ). X ( { } 1 2 , ,..., X N X X X X = with X N
The model Z can be replaced by an approximate model (metamodel) ( )

ˆ, a LAR Z f X =
and is given in the following equation:

( ) ( ) , a , a . LAR Z f X Z f X = ≈ = (5)
The metamodel will be discussed in the next section. The Monte Carlo Simulation is applied to the model ( ) Z X or ˆ( ) Z X in order to estimate the POD.

• Stochastic Metamodel

It is necessary to find a model faster than the FV model in order to reduce the calculation time. In this case, it may be considered to propagate the uncertainties through a metamodel constructed from Truncated Polynomial Chaos Expansion [START_REF] Ghanem | Stochastic Finite Elements. A Spectral Approach[END_REF]. In that case, the approximation ˆ( ) Z X of the stochastic FV model ( ) Z X can be written in the form [7]:

( ) ( ) 1 ˆ( ) , out P i i i Z X X α ξ = = Ψ ∑ (6) N 3 N 1 N 2 N 4 P n r c ∆ x y η ∆
) ! ! ! X out X N p P N p + =
, if the number of the input parameters is 6

X N = and the polynomial order 3 p = , thus the number of polynomial is equal to 84 out P = . The coefficients of the polynomial can be estimated using a non-intrusive method such as the regression approximation [START_REF] Lemaître | Spectral Methods for Uncertainty Quantification[END_REF]. The number of realisations desired is at least out P , therefore, X out N P ≥

. The unknown coefficients i α can be calculated by the least square minimization, i.e. by minimizing the mean square of the residual reads:

( )

2 1 ârg min E ( ) ( ) , X N j j j Z X Z X α =     = -     ∑ (7)
with the mean operator [ ] E . . The LAR method is a regression method that reduces the computational cost by selecting the polynomial terms which are correlated the most with the output. The number of terms is therefore significantly reduced compared to the classical regression method. However, the number of polynomials increases exponentially with the dimension. So, in order to decrease the number, an improved LAR which consists in constructing iteratively with an increasing order of polynomials has been used [START_REF] Nguyen | Uncertainty Quantification Using Sparse Approximation for Models With a High Number of Parameters, Application to a Magnetoelectric Sensor[END_REF].

• Threshold Determination

To determine the detection limit (s), the probability of False Alarm (PFA) has to be imposed at a very low value [START_REF] Rajesh | Numerical Model Based Approach for Estimating Probability of Detection in NDE Applications[END_REF]. In our case, we put PFA = 0.05. Therefore the detection limit s is calculated by the following formula: 0 ( ) 0.05.

PFA Z s = > = P ( 8
)
with 0 Z is the random impedance without any defect.

• Probability of detection

After the detection limit has been determined, the POD is estimated. The POD is the probability of the model response Z of a given defect ( a 0 > ) that is above threshold s as shown in the following expression

( ) ( ). POD Z s Z s = > ≈ > P P (9)
The total number of realizations of the impedance in the presence of defect is denoted by M , we denote s M the number of realizations that exceeds the threshold s , therefore the POD can be estimated by / s POD M M = .

IV. GEOMETRIC TRANSFORMATION (GT) METHOD

The dimensions of the sensor as well as the thickness of the tube are random input parameters. It means that during the sampling process the geometry is modified. Due to the large number of FV model calls, the remeshing of the geometry adds non negligible computational time. Moreover, the remeshing introduced an additional numerical noise due to the modification of the connectivities between nodes. The GT method is based on changing the position of the coordinates of the nodes without changing their connectivity. First of all an initial mesh is proposed and decomposed into different subdomains, if the coordinates are changed according to the random vector g X (vector of geometric random parameters with g X X ⊂ ), then the new mesh is obtained by dilation, compression or translation. The main idea of the GT method is to determine the appropriate transformation. Fig. 3 represents an initial domain and a transformed domain knowing that the new coordinates of g X .

The transformation T transforms the initial domain E into a random domain ( )

g D X .
For each realization of X, the nodes are repositioned. Fig. 4 gives the different steps to obtain with the two methods (GT and RM) the stochastic approximate model of the impedance by applying the LAR method.

V. NUMERICAL APPLICATION

A sample of M realizations of the input parameters is generated. With the RM method the geometry and the mesh are recalculated for each realization. Therefore, we estimate the POD for different meshes with a different size, denoted E ( ) The second method used is the GT method, in which the mesh is fixed at the beginning (initial mesh) and the appropriate transformation is determined. Meshes of different size are also considered as in RM method (S1, S2, S3, and S4). The process of estimating the POD for the GT method denoted POD GT . For each realization, the GT is used in order to relocate the nodes according to the modification of the geometry.

g D X ( , ) g T T x X = (b) (a)
In the both methods (RM and GT), the POD RM and POD GT will be determined with defect sizes a 1 , a 2 , and a 3 equal to 0.05, 0.4970, and 0.9675mm respectively.

To estimate the POD with two methods, the approximate model is used. It is constructed with 150 N = realisations of the FV model obtained with GT or RM method, once the metamodel is constructed, the MCSM is applied with a sample of 10000 M = realisations for each value of defect size a i . Repeating the same process with the RM and GT methods and all meshes. A linear metamodel is constructed and considered in this study; hence, the chaos coefficients are perfectly calculated.

To verify the accuracy of the approximate model, an error is calculated by the following expression:

( ) ( ) 1 1 M j j j Z X Z X ε =     = -     ∑ (10)
The convergence of the metamodel is accurate with an error less than 10 -4 which shows that the approximation is of good quality.

The Fig. 5 ; it shows great agreement with the two methods (RM and GT) based on LAR and the reference.

The Table I represents the POD for different methods (RM and GT) and for different size meshes (S i ) where we can see also a good agreement between the different methods with changing of the defect size a i . We can see that accurate results can be obtained with a coarse mesh and we can see also that the RM and GT methods give very similar results for this case meaning that the numerical has really few influence in the considered example.

VI. CONCLUSION

In this paper two approaches to estimate the POD of an eddy current testing problem of tube coil system (steam generator tubes) have been estimated. The modification of the geometry which is random has been handled by a geometric transformation. The geometric transformation (or mapping) is shown to give similar results to the remeshing technic but is less time consuming. The stochastic metamodel has been used in both methods to reduce the computation cost of the high number of realisation and showed a great precision.

Fig. 2

 2 Fig. 2 Unknown elements representation.

(

Fig. 3 Fig. 4

 34 Fig. 3 Geometric domain; a) initial domain ( E ); b) transformed domain ( ( ) g D X ).

  represents the POD for the S2 mesh, the POD Ref is the reference POD of the GT method with 10000

Fig. 5

 5 Fig. 5 Probability of Detection for different methods (RM and GT).


Jean-François Katzensteiner
email: klaus.katzensteiner@boku.ac.at

B Ponge

G Jabiol

E Sartori

J.-F Ponge

Le Bayon

TerrHum: An iOS App for Classifying Terrestrial Humipedons and some Considerations about Soil Classification

  L'archive ouverte pluridisciplinaire

Introduction

In 2003, 26 soil scientists gathered in Trento (Italy) to standardize methods and rules for classifying the biologically most active upper part of the soil. For doing so, it was essential to develop a common terminology and to exchange data on soil carbon dynamics and its relation to morphological characteristics of organic and organic-mineral topsoil. The definitions of organic and organic-mineral horizons of the soil, for example, diverged a lot. These scholars came from 9 different European countries and immediately realized that the task was not an easy one because each country had its own historical classification systems, often incompatible in-between. It took more than 5 years of meetings and discussions to publish, on behalf of a part of the members of the 'European Humus Group", then simply 'Humus Group", a first synthesis [START_REF] Zanella | Toward a European humus forms reference base[END_REF].

The main problem was not the existing classifications, but the fact that for bringing together the existing pieces of knowledge to form a frame compatible with all points of view, it was necessary to enlarge the picture, and also to understand what the different classifications had not considered individually. As an example, there were some humus forms called 'humus gemellare" (twin humus) by an Italian classification. In France, these forms were considered as 'double" forms, but rather classified as Mulls (Amphimull) because of their thick and impressive crumby organic-mineral horizon [START_REF] Jabiol | Une clasification morphologique et fonctionnelle des formes d'humus. Proposition du Référentiel Pédologique[END_REF]. In other countries, these same forms were classified as Moder, because of their very thick OH horizon [START_REF] Broll | Topsoil characterization -Recommendations for revision and expansion of the FAO-Draft (1998) with emphasis on humus forms and biological features[END_REF].

Because of all these disparities, it was then proposed to bring all these 'double" forms together into a new category called 'Amphi". This decision caused a rift between i) those who had seen and supported the existence of such forms of humus, and ii) those who instead preferred to stick to the current classification systems without considering such forms, because they were not very common in their respective countries. The issue was solved only after a workshop in San Vito (2005, Italy), with field trips to see 'Amphi" and comment on it.

Back home, a German member of the Humus Group published an article explaining to his compatriots the existence of such forms in Germany, also highlighting that they could become more frequent and therefore be used as indicators of ongoing global warming [START_REF] Graefe | Gibt es in Deutschland die Humusform Amphi? Mitteilungen der[END_REF].

However, the story went on. When members of the Humus Group came together to see these Amphi forms in Austria (Vienna, 2004), other humus forms already described many years before did not match this category, and then fell more or less in oblivion [START_REF] Hartmann | Forstökologie. Zustandserfassung und Standortsgemässe Gestataltung der Lebengrundlagen des Waldes[END_REF][START_REF] Hartmann | Waldhumusdiagnose auf Biomorphologischer Grundlage[END_REF][START_REF] Kubiëna | The Soils of Europe. Illustrated Diagnosis and Sistematics[END_REF]. So are Tangel forms (sometimes called 'Alpenmoder" by local forest managers), found at high altitude on calcareous substrates, while corresponding Mor humus forms, with different morphological and chemical properties, are related to siclicate bedrock and acidic soils. The fact opened a diatribe concerning differences between Mor and Tangel that lasted for years and was solved only recently [START_REF] Kolb | Tangel humus forms -genesis and co-evolution with vegetation[END_REF][START_REF] Meynier | Humus forms pathways in low-elevation cold scree slopes: Tangel or Mor?[END_REF]. In order to clarify the situation, we decided to set up a classification linked to the soil forming substrate (only if it influences the topsoil: Amphi and Tangel on calcareous substrates, Moder and Mor on 'acidic' ones, with Mull in-between in neutral or nearly neutral situations. All members of the Humus Group were then in agreement with this solution, except supporters of 'Mor rather without pedofauna" (French point of view) and others who wanted instead 'Mor with rather relevant pedofauna" (German supporters who described Mor as extreme Moder forms).

One day, inevitably, the classification bogged down on the definition of soil (Zanella et al., 2018a;e;[START_REF] Zanella | Editorial. Humusica 1 -Terrestrial Natural Humipedons[END_REF] a point of contention for more than two years. For some of us, abiotic factors predominate pedogenesis, while for others biotic factors outshine.

When the soil was proposed as a 'digestion system", a 'belly", with living organisms having it inside (in their belly), the group started to crumble. And it spread to individual members when the idea came out that soil binds life and death (for many soil scientists, philosophy should remain separated from what they call 'real science", 'numbers", 'mathematical evidence", even if all of us -as living beings -are destined to die an incalculable day coming). Against all odds, however, the group remained united, but encountered very often many other difficulties, for example when the various humus forms were arranged on a global scale (Zanella et al., 2018f;[START_REF] Ponomarenko | A recommendation for a new descriptor for pyrogenic soil horizons in the Canadian Soil Classification System[END_REF] . It was then necessary to introduce the concept of 'humus system", comprising ecologically and functionally similar humus forms as well as creating new systems for collecting strange, unusual or never surveyed humus forms (Zanella et al., 2018d) or for those in water (Zanella et al., 2018l, k;b) and in anthropogenically strongly shaped environments (Zanella et al., 2018j).

Nevertheless, there are still many questions that may be resolved in the near future. New methods at the frontier of science in soil biology can answer issues linked to co-evolution of soil and organisms, biogeography of soil organisms, soil biodiversity etc.. Anyway, all these questions give a huge importance to the soil and may be reflected in its morphological characteristics. Humans need to better know the soil, but what about 'humus'? Recently it has been suggested to abandon the use of this term, because 'humus' cannot be clearly characterized from a chemical point of view [START_REF] Lehmann | The contentious nature of soil organic matter[END_REF]. We strongly suggest to keep the term and to speak of 'humus' when organic particles and accompanying biota are present.

How does TerrHum work?

Let us consider a user facing a soil profile to be classified. A cubic 50x50x50 cm hole dug into the ground is generally sufficient for studying humus systems and forms in a forest environment, while a larger hole or many well-distributed small holes are necessary for a representative survey of a heterogeneous area (Zanella et al., 2018i). A detailed description of the main actors of soil biodegradation and their relationships with the horizontal and vertical distribution of humus horizons is reported in (Zanella et al., 2018c).

As a humus form is made of superposed humus horizons, the app asks the user to indicate one by one which types of humus horizons are present in the observed profile. Organic horizons (OL, OF and OH, corresponding to USDA Oi, Oe, Oa) were clearly distinguished from organic-mineral ones (different types of A horizons) (Fig 1).

In general, when in the field a horizon appears totally consisting of organic remains, it is classified as an organic horizon. It is well accepted that the organic matter (OM) composing such horizons amounts to more than 1/3 of the total weight, which corresponds to 20 % or more organic carbon (OC). That allows the separation of organic and organic-mineral horizons in the field and the adjustment of the classification, if necessary, after laboratory analyses of the OC concentration in a horizon sample. The organic horizons are divided into OL, OF and OH following their decreasing content of "recognizable remains" (e.g. a leaf, a needle, a piece of bark…): more than 90 % in OL, from 90 to 30% in OF, and less than 30% in OH.

Another diagnostic criterion is, if the process of litter transformation (decomposition) is accomplished mainly by animals or by microbes. Respectively, the process generates zoogenic zoOF or non-zoogenic nozOF horizons described in the app and recognizable in the field.

Accordingly to its structure and genesis, the organic-mineral A horizon is classified in zoogenic biomacro, biomeso and biomicro A horizons (maA, meA, miA), and in non-zoogenic massive and single-grain A horizons (msA and sgA). All these horizons are described and illustrated in the app.

The user is asked to answer a series of YES/NO questions in a dichotomous key; an example: 'is OH horizon present?" (Fig. 2,3a,3b,3c,3d). Slightly differing from the App, the simplified key (Fig. 2) requires some field experience, but allows a faster, equivalent, correct classification. For detailed definitions of all diagnostic horizons and criteria of classification, refer to Zanella et al. (2018l, 2018m).

A touch-button located at the bottom of the screen allows the user to recall at any time definition and photographs of each diagnostic horizon, allowing users to more accurately define the real horizon (Figs 4a,4b,4c,4d).

Other keys at the bottom of the screen allow retrieving examples of humus systems and forms (Fig. 5), as well as tables of composition/classification of humus horizons, or groups of animals and their droppings (Fig. 6). At the end of the classification process, a photograph of the target humus form appears, along with a list of the chosen horizons.

By touching the screen, each photograph may be magnified (Figs 3d,4b,5b,6b). A caption at the bottom of each picture provides access to the morpho-functional features of the soil profile and leads finally to the classification of the humus system.

Many examples of soil horizon features are documented in the app. Future updates of the app will include useful links to enriched external data banks. We are currently collecting information and photographs to further improve the classification of charcoal, charred mass and other dark-coloured sublayers that can interfere with the definitions of standard diagnostic horizons (e.g. [START_REF] Ponomarenko | A recommendation for a new descriptor for pyrogenic soil horizons in the Canadian Soil Classification System[END_REF]). We also are collecting new photographs and drawings of soil animals, related to the assessment of soil biological quality.

Is it possible to integrate TerrHum into different soil classification systems?

The purpose of TerrHum is to share a morpho-functional classification of topsoils at global level. All over the world, students, researchers, forest managers may recognize humus forms, take a georeferenced photograph and send it to one of the authors of this article. Collected and checked data may help to improve published maps (Zanella et al., 2018f) or to prepare new soil and humus maps. The photographs being georeferenced, it is relatively easy to create new maps if the number of points is sufficient.

Fire is an integral process of carbon transformation and pyrogenic features of humus horizons are currently overlooked, and still have to be incorporated into the TerrHum system. Adding these features to the Terrhum app may enable collecting an important layer of information on past fire events at both local and global levels.

An Android version of TerrHum for smartphones will be available in the near future.

TerrHum (iOS) and/or the up-coming Android version or in general the concept of TerrHum can be an example also for soil classification. For this purpose, it would be important to adopt a soil model subdivided into Humipedon, Copedon and Lithopedon published in (Zanella et al., 2018g;[START_REF] Zanella | f. Humusica 2, article 18: Techno humus systems and global change -Greenhouse effect, soil and agriculture[END_REF] We programmed to work on apps like this for Copedons and Lithopedons. Histosols are also considered as humus forms and will be later integrated in TerrHum. People interested in the project are kindly asked to contact the authors of this paper. Given Humipedons, Copedons and Lithopedons could then be associated to provide a full profile and name for each type of soil. This is the reason we would like to widen the concept of soil. We would like to divide the soil profile into three parts, each of them being studied either separately or brought together to better understand the overall soil functioning (Figs. 7 and8). We know that this concept may hurt some among the pedologists, however, it is crucial to understand that each humus form functions relative to soil and vegetation but has also its own spatio-temporal scale of formation, functioning and dynamics [START_REF] Bernier | Humus form dynamics during the sylvogenetic cycle in a mountain spruce forest[END_REF].

The soil may be divided into three layers which may be described and classified relatively independent each from each other: Humipedon, Copedon and Lithopedon. The first depends mainly on the source of organic matter and animals and microorganisms that live in the soil: they generate the Humipedon composed of organic (OL, OF and OH) and organic-mineral (A) horizons [START_REF] Hole | Effects of animals on soil[END_REF].

Beyond these, in the most evolved soils, different mineral horizons like E and B can be distinguished. Their development, though linked to root dynamics and turnover of organic matter is strongly dependent on physical and chemical soil processes. As for those belonging to the Humipedon, there are different types of Copedon horizons, according to climate and parent material [START_REF] Muhs | Impact of Climate and Parent Material on Chemical Weathering in Loess-derived Soils of the Mississippi River Valley[END_REF].

The bottom of the profile completely depends on the original rock. These are the fragmented C horizons and the hard R rock layer. There are different types of Lithopedon (e.g. carbonatic, silicatic; all possible morpho-functional types or more or less transformed initial rocky basements). The Lithopedon is not always at the bottom of a well-differentiated soil profile.

A very young soil is just a thin 'biofilm" laying on a Lithopedon, like in deglaciated environments [START_REF] Wynn-Williams | Response of pioneer soil microalgal colonists to environmental change in Antarctica[END_REF]. Then Humipedon and Lithopedon are built-up, the soil profile evolving toward a complete sequence of Humipedon, Copedon and Lithopedon (Fig. 5). Occasionally the soil can either lose some of its more superficial parts by erosion (truncated soils; [START_REF] Desmet | GIS-based simulation of erosion and deposition patterns in an agricultural landscape: a comparison of model results with soil map information[END_REF]) or one part grows till taking all the place (e.g. Humipedon in submerged soils -Histosols, or Copedon in tropical Vertisols).

Once mastered this classification, things become simpler and easier to understand; and the definitive/detailed name of the soil could arise from combining the names given separately to the humipedon and the classification of the soil profile in any system. The subdivision into strata of the soil profile. The biodiversity maximun is not realized at the surface, but in the layer that separates the horizons OL and OF or in the horizon OF. The latter is a very protected horizon (in particular desiccation but also rain), which contains a large number of trophic resources (very transformed dead leaves, but also fungi and bacteria) and which, by its fragmentation, constitutes many refuges for species that would otherwise be in strong competition, all that is needed to have maximum biodiversity. Figure 7. The subdivision into strata of the soil profile. The biodiversity maximun is not realized at the surface, but in the layer that separates the horizons OL and OF or in the horizon OF. The latter is a very protected horizon (in particular desiccation but also rain), which contains a large number of trophic resources (very transformed dead leaves, but also fungi and bacteria) and which, by its fragmentation, constitutes many refuges for species that would otherwise be in strong competition, all that is needed to have maximum biodiversity.

List of Figures:

Figure 1 .

 1 Figure 1. The basic concept of humus classification and the designation of organic layers on

Figure 2 .

 2 Figure 2. Dichotomous key for identifying Terrestrial Humus systems and Forms. The

Figure 3 .

 3 Figure 3. TerrHum screens. a) Initial screen. The main key of classification is accessible by

Figure 4 .

 4 Figure 4. TerrHum screens. a) Clicking on the first touch-button 'Horizons" (Fig. 3a) allows to

Figure 5 .

 5 Figure 5. TerrHum screens. a) Clicking on the third touch-button'Systems and Forms" (Fig. 1a)

Figure 6 .

 6 Figure 6. TerrHum screens. a) Selecting the fourth touch-button 'Tables" (Fig. 3a) opens the

Figure 7 .

 7 Figure 7. The subdivision into strata of the soil profile. The biodiversity maximun is not realized

Figure 8 .

 8 Figure 8. Pedon sub-units: Humipedon, Copedon and Lithopedon in a theoretical series of soil

Figure 1 .

 1 Figure 1. The basic concept of humus classification and the designation of organic layers on forest floors (O-horizons) and organic-mineral soil horizons in the European classification of humus systems (TerrHum). 338x190mm (96 x 96 DPI)

Figure 8 .

 8 Figure 8. Pedon sub-units: Humipedon, Copedon and Lithopedon in a theoretical series of soil development. From (Zanella et al., 2018e; a), slightly modified. 81x60mm (200 x 200 DPI)

Soil Sci. Soc. Am. J. Accepted Paper, posted 03/21/2019. doi:10.2136/sssaj2018.07.0279