Raphaël Prat
email: raphael.prat@cea.fr

Laurent Colombet
email: laurent.colombet@cea.fr

Raymond Namyst
email: raymond.namyst@u-bordeaux.fr

Combining Task-based Parallelism and Adaptive Mesh Refinement Techniques in Molecular Dynamics Simulations

Keywords: Adaptive Mesh Refinement, Task Parallelism, Molecular Dynamics

Modern parallel architectures require applications to generate massive parallelism so as to feed their large number of cores and their wide vector units. We revisit the extensively studied classical Molecular Dynamics N-body problem in the light of these hardware constraints. We use Adaptive Mesh Refinement techniques to store particles in memory, and to optimize the force computation loop using multi-threading and vectorization-friendly data structures. Our design is guided by the need for load balancing and adaptivity raised by highly dynamic particle sets, as typically observed in simulations of strong shocks resulting in material micro-jetting. We analyze performance results on several simulation scenarios, over nodes equipped by Intel Xeon Phi Knights Landing (KNL) or Intel Xeon Skylake (SKL) processors. Performance obtained with our OpenMP implementation outperforms state-of-the-art implementations (LAMMPS) on both steady and micro-jetting particles simulations. In the latter case, our implementation is 4.7 times faster on KNL, and 2 times faster on SKL.

INTRODUCTION

The High Performance Computing community is currently witnessing a significant increase of the number of cores in parallel supercomputers, which are already featuring more than 10 millions of cores. This increase of the parallelism degree is a clear trend that unfortunately comes with a decrease of the amount of memory per core, as well as a strong pressure on the use of highly vectorized code. The Intel Knights Landing processor (knl) [START_REF] Sodani | Knights landing (KNL): 2nd Generation Intel® Xeon Phi processor[END_REF] and Intel Skylake processor (skl) [START_REF] Doweck | Inside 6th-Generation Intel Core: New Microarchitecture Code-Named Skylake[END_REF], for instance, require respectively from 64 to 256 and 48 to 96 threads to utilize the cores at their full capacities, as well as 512-bit wide SIMD instructions (i.e. 16 floats). This explains why many ongoing efforts are devoted to designing new implementations of parallel simulations to meet these hardware constraints, which require to generate more parallelism, more vectorized code and to enforce stronger data locality.

Molecular Dynamics (MD) N-body simulations, which represent the most accurate method to simulate complex materials by operating at a microscopic level, are naturally good candidates for exploiting most of the underlying hardware's potential, because they are intrinsically data-parallel. Optimization of such simulations over parallel supercomputers has indeed been extensively studied, and very efficient implementations are available over a wide range of modern hardware. One of the key aspect of an efficient implementation of an N-body simulation is to speedup, for each particle, the lookup of neighbor particles which potentially exerts a force on the considered particle. To this end, several methods are used by state-of-the-art software, most notably the Verlet list [START_REF] Verlet | Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF] method or the cell-list one [START_REF] Allen | Computer Simulation of Liquids[END_REF]. These methods perform indeed very efficiently under regular conditions, that is, when particles are uniformly spread over the domain. For extreme shock simulations where many particles have high velocities, these methods lead to either frequently recompute neighbor lists or to cope with a massive number of empty cells.

This paper investigates the use of Adaptive Mesh Refinement techniques [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] in N-body simulations. In order to cope with regions of highly variable particle density, we partition the domain in several subdomains organized as octrees, which are periodically refreshed. Each octree is intended to be processed sequentially on a single core, and we carefully enforce the execution order so that no adjacent octrees can concurrently access a shared cell, leading to a completely lock-free algorithm. To maximize the efficiency of particles interactions computation over recent manycore architectures, such as Intel Xeon Phi knl, we use a careful data layout for particle attributes that improves data locality and cache utilization, and maximizes code vectorization opportunities. We evaluate our approach using an extended implementation of the ExaStamp MD simulator [START_REF] Colombet | ExaStamp: A Parallel Framework for Molecular Dynamics on Heterogeneous Clusters[END_REF] over a Intel Xeon Phi and Intel Xeon. We use different input configurations and show that our approach outperforms state-of-the-art solutions while using a similar memory footprint.

The main contributions of this paper are:

• We introduce a new software architecture for MD including amr techniques; • We present a thread parallelization strategy adapted to the amr structure; • We evaluate our implementation of the amr on two different architectures using both a static and a highly dynamic test case.

The remainder of this paper is organized as follows. In the next Section, we briefly recall the main concepts used in Molecular Dynamics simulations and in Adaptive Mesh Refinement based applications. Section 3 describes the main ideas behind our approach.

We demonstrate the relevance of our approach in Section 4 comparing our implementation with state-of-the-art software on several benchmarking scenarios. We present related work in Section 5, and dress concluding remarks in Section 6.

BACKGROUND 2.1 Molecular Dynamics Simulations

The classical MD is a computational method used to study dynamical properties of a systems made of particles. This method consists in numerically integrating Newton's equation of motion

- → f = m. - → a 1 ,
where the force on a particle depends on the interactions with all others [START_REF] Allen | Computer Simulation of Liquids[END_REF]. The force on particles is given by the gradient of a potential

- → f =-∇ - → U.
MD simulations typically follow an iterative process, where each iteration (time step) consists in computing forces applied to each particle, deducing each particle's acceleration, updating particles velocities, and finally updating particles positions. To do so, the most widely adopted scheme is the Velocity Verlet integrator [START_REF] Benedict | Integration Methods for Molecular Dynamics[END_REF].

Particle interactions.

In most MD simulations, particles are processed as points and the interacting force between particles is approximated as a gradient of a potential that depends on the particles relative positions. The force computation is obviously the most challenging part: it contains all the physics of the simulation and can take up to 95% of the total time. Potentials from classical MD are empirical or semi-empirical, computed from an analytical formula or interpolated from tabulated values [START_REF] Wolff | Tabulated Potentials in Molecular Dynamics Simulations[END_REF].

In this paper, we focus on short-range interactions neglected beyond a given distance rcut called the cutoff distance. This approximation is relevant for solid materials, in which distant atoms are "screened" by closer atoms.

Lists of Neighbors.

To speed up the force computation step, most MD simulations maintain a list of neighbors for each atom. To build such lists efficiently, two methods (see figure 1) are typically used and can even be combined. The first is based on Verlet lists [START_REF] Verlet | Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF]. A radius of Verlet (rVerlet) is added to rcut. As a consequence, the list of neighbors contains "useless" atoms which are too far from the considered atom. However, as long as no atom has moved from its original position 2 by more than 1 2 rVerlet, 2 The original position is the one captured when updating the list of neighbors. there is not need to recompute the lists. This method is therefore the most efficient when atoms are moving very slowly, because the Verlet lists are not frequently updated. This method have a complexity of O(N 2 ), and can easily be performed in parallel. The other one is the linked cells method [START_REF] Allen | Computer Simulation of Liquids[END_REF]. The domain is divided into rcut sized cells on a grid. For all atoms, the neighborhood of an atom is included in the (27 in 3D) neighboring cells. This method avoids the cost of maintaining lists of neighbors per atom. However, it comes at the cost of considering a significant amount of atoms which are beyond the rcut distance. This method has a complexity of O(N ). The strategy of building neighbor lists can influence the simulation performance because it requires a lot of memory accesses.

1 - → f force, m mass, - → a acceleration, - → U potential energy

Pair Potentials.

Pair potentials are a particular but extremely widely used case of potentials, where interactions between particles are reduced to pairwise interactions. Therefore, the potential between two particles will only depend on the distance between those particles. Pair potentials are quite effective to describe lowdensity materials such as liquids and gases. Some common pair potentials are for instance, the Lennard-Jones potential (lj) [START_REF] Jones | On the Determination of Molecular Fields. II. From the Equation of State of a Gas[END_REF], the Morse potential [START_REF] Morse | Diatomic Molecules According to the Wave Mechanics. II. Vibrational Levels[END_REF] and the Exponential-six potential [START_REF] Mason | Transport Properties of Gases Obeying a Modified Buckingham (Exp-Six) Potential[END_REF]. Although it was first designed to study gases, the Lennard-Jones potential (lj) has been used in a large part of materials science and has become a standard benchmark for MD codes.

Adaptive Mesh Refinement (AMR)

The first adaptive mesh refinement was developed by Berger [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] for hydrodynamics application in 1984 to optimize time computation. In hydrodynamics, the domain is discretized into a Cartesian grid where the solution is computed in each point of the grid. The spacing of the grid points determines the local error and hence the accuracy of the solution. amr consists in coarsening the grid points where the points do not need to be spaced to keep their stability. The error is negligible at these points. The coarsening is determined by a criterion based on density or on a Richardson extrapolation [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Fry | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF]. Note that different levels of refinement may be obtained for a simulation. This method is interesting for simulations with heterogeneous density to reduce computational cost for some parts of the simulation. amr methods have already been introduced in N-body simulations in the past. In such simulations, the domain is decomposed into a grid of cells, which are equivalent to points in hydrodynamics. The amr consists in coarsening the cells containing only a few atoms, using criteria such as the number of particles or average density. A large amount of amr software have been developed [START_REF] Dubey | A survey of high level frameworks in block-structured adaptive mesh refinement packages[END_REF] for various physical problems such as astrophysics with Enzo [START_REF] Greg L Bryan | Enzo: An adaptive mesh refinement code for astrophysics[END_REF] or BoxLib [START_REF] Lijewski | Boxlib[END_REF], Chombo [START_REF] Colella | Chombo software package for amr applications-design document[END_REF] to solve PDE.

ADAPTIVE MESH REFINEMENT FOR MOLECULAR DYNAMICS SIMULATIONS

To meet the needs of highly dynamic particles simulations over many core architectures, we introduce a new approach inherited from Adaptive Mesh Refinement Methods (amr) to perform MD simulations. We first describe the main characteristics of the Ex-aStamp framework (Section 3.1) that has served as a basis for our implementation. Section 3.2 presents the main data structures used to efficiently cope with highly heterogeneous sets of particles. In Section 3.3, we detail how the workload is dynamically assigned to the underlying computing resources.

The ExaStamp MD Framework

ExaStamp is a MD code which has been recently developed at CEA for the upcoming generation of supercomputers. The main goal of ExaStamp is to cope with simulations featuring up to several billions of atoms on thousands of cores. The code is written in C++ and uses hybrid parallelization mixing MPI and threads (Open-MP/ Intel TBB). ExaStamp kernels are written by physicists using high-level types and functions, and the generated code features Intrinsics instructions to minimize missed vectorization opportunities. These techniques have most notably proved to be effective on the Intel Xeon Phi Knights Corner processor.

Like many MD software packages, ExaStamp relies on the linkedcell method to organize particles, that is, a static mesh of cells that does not evolve during the simulation. Although this strategy exhibits high performance when particles are homogeneously distributed across the simulation domain, it is not suitable for extreme dynamic micro-jetting [START_REF] Durand | Comparative simulations of microjetting using atomistic and continuous approaches in the presence of viscosity and surface tension[END_REF] configurations. Indeed, between two time steps, most atoms move from one cell to another, leaving a large number of cells almost empty. Dealing with those empty cells incurs significant memory and performance overheads, because their number can sometimes exceed the number of atoms.

Software Architecture

Our approach combines the use of several techniques such as the octree's structure, the Structure-of-Array, the Morton index and a vector-friendly handling in order to optimize access to L1 and L2 caches.

Adaptive Mesh

Refinement. We apply amr techniques to the management of cells, so as to deal with cells of different sizes, depending on the particles density. Figure 2 illustrates how the spatial domain is decomposed into cells [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. The domain is first decomposed into sub-domains, each sub-domain being assigned to one MPI process. Sub-domains are in turn decomposed into a 2) is decomposed into sub-domains that are distributed across MPI nodes, each sub-domain being a grid of octrees. grid of octrees. The reason for using a grid of octree instead of a large, unique octree is that an octree serves a unit of work that can be exchanged between MPI processes for load balancing purpose. octrees are built and refreshed using a refinement algorithm: each cell is recursively divided into 8 cells of same volume if some refinement criteria are fulfilled. The maximal depth D max limits the number of times octrees can be refined. Therefore each octree contains between 1 and 8 D max cells. As shown in Figure 2, a cell can either be a root cell (also called level 0 cell), a leaf cell or an inter cell. Note that the minimal size of a leaf cell corresponds to the size that would be used in a classical linked-cell method, in accordance with the maximum cut-off distance.

The way octrees are assigned to MPI processes is beyond the scope of this paper. Inside MPI processes, thread-level parallelism is used to explore the grid of octrees, as it will be detailed in the following Subsection. Given that the number of octrees is larger than the number of cores by several orders of magnitude, one important aspect of our strategy is that each octree is assigned to a single thread and is thus processed sequentially. This allows us to optimize the octree traversal with respect to cache usage [START_REF] Gannon | Strategies for cache and local memory management by global program transformation[END_REF][START_REF] Dennis | The influence of memory hierarchy on algorithm organization: Programming FFTs on a vector multiprocessor[END_REF]: leaf cells are explored along a Z-order curve to minimize cache misses, as illustrated in Figure 3 (red path).

Data layout.

Inside an octree, all atoms' attributes are stored in the root cell, in a Structure-of-Arrays (SOA) layout. For the sake of clarity, we simply refer to this data structure as the "array of atoms" in the remaining of the paper. inter and leaf cells only store indexes indicating the position of atoms in this array, as illustrated in Figure 4. This strategy prevents the copy of atoms during the refinement, improves the data locality between cells inside an octree and reduces data movements because transfers are performed between octrees instead of cells [START_REF] Durand | A packed memory array to keep moving particles sorted[END_REF]. In addition, because the number of arrays is reduced, they are filled with more atoms, which improves vectorization efficiency.

For each atom, we compute a Morton index [START_REF] Guy M Morton | A computer oriented geodetic data base and a new technique in file sequencing[END_REF] out of the cartesian position of its smallest bounding cell. This index gives the path from the root cell to the leaf containing the atom. The array of atoms is sorted by the atoms' Morton index, to improve cache locality [START_REF] Meloni | Efficient particle labeling in atomistic simulations[END_REF]. This strategy avoids the copy of data during refinement/coarsening steps.

Vector-friendly

Particles Handling. To maximize vectorization opportunities inside the computation step, one technique consists in using a vectorization buffer to pack information about neighboring atoms in contiguous areas of memory. While this allows the compiler to generate efficient, vectorized instructions to process data stored in such buffers, the cost of introducing an extra copy of data has a negative impact on performance in the case of low cost potentials (e.g. Lennard Jones). Another technique is to access the data through in a vector-friendly manner. More precisely, we divide the set of atoms into a series of blocks (for which the size is a multiple of the processor's vector width), which is illustrated in Figure 5. We only process blocks that contain at least one atom included in the Verlet radius, while others are skipped. In the worst case, we could end up with blocks containing only one useful atom, which would lead to process only one single atom at a time. In practice, however, we observe that vector registers are filled with more than 50% of useful values. This implementation variant of ExaStamp will be called amr Vec in the remaining of this paper.

Parallelization Strategy

Inside each MPI node, we use OpenMP to process cells in parallel. As mentioned previously, since the number of octrees is significantly larger than the number of cores, each octree is assigned to a single OpenMP thread, and thus is processed sequentially. In our current implementation, we spawn one OpenMP task per octree.

When the potential has symmetrical properties (third newton's law), forces exerted between two particles are equal and are thus computed once only. This saves some precious computing cycles, but requires as a counterpart to use costly synchronization tools (e.g. mutexes) to avoid that multiple threads update the same target particle. Another strategy is to process the cells in an order that ensures that no particle can be simultaneously updated by different threads. The wave method was introduced by Meyer et al. [START_REF] Chris | A hybrid algorithm for parallel molecular dynamics simulations[END_REF][START_REF] Meyer | Efficient parallelization of molecular dynamics simulations with short-ranged forces[END_REF][START_REF] Meyer | Efficient parallelization of short-range molecular dynamics simulation on many-core system[END_REF] as a mean to avoid these synchronizations. The idea is to partition cells -or octrees in our case -into a set of waves such that two cells of the same wave are far enough from each other to have no neighbor in common. Thus, during the computation of a given wave, neighbor particles can be updated without synchronization (see Figure 6).

Note that there is no need to strictly execute waves one after the other: their execution can actually overlap if it relies on a task-based scheduler able to take inter-task dependencies into account (such as Intel TBB or OpenMP). During the execution of wave n, tasks of wave n + 1 will be progressively unlocked. Using dependency graph allows to unlock tasks as quickly as possible. We use OpenMP 4.5 to generate explicit tasks with dependencies. In our case, each task represents an octree to compute (see Figure 7), and not just a cell. This method is equivalent to the cell task blocking method [START_REF] Meyer | Efficient parallelization of short-range molecular dynamics simulation on many-core system[END_REF] and limits the number of waves to 8 but requires that D max is superior to one. This approach is similar to the Particle and Separate-Calculation (PSC) [START_REF] Hu | Kernel optimization for short-range molecular dynamics[END_REF].

Finally, compared to the original version of the wave method, octrees helps to control the granularity of tasks, and thus helps to minimize the overhead of generating OpenMP tasks.

EVALUATION

To assess the effectiveness of our approach, we present a performance study using two test cases: bulk material, a steady configuration of atoms, and micro-jetting, a highly dynamic one. We compare our implementation (ExaStamp extended with amr) against the Lammps [START_REF]LAMMPS Molecular Dynamics Simulator[END_REF] reference package and its mini-application min-iMD [START_REF] Michael | Improving performance via mini-applications[END_REF].

Experimental Framework

Our experimental framework is composed of two architectures. The first is a knl processor featuring 68 cores running at 1.4GHz (4 of which being dedicated to the operating system). Each core has a 2-wide out-of-order pipeline with support for 4 hardware threads and 512-bit vector units. Its memory system consists of a 96 GB DDR4 off-chip memory and a 16 GB MCDRAM stacked memory. The second architecture is a skl-based machine that contains 2 sockets of 24 cores running at 2.7GHz. Each core has a 2-wide out-of-order pipeline with support for 2 hardware threads and 512-bit vector units. It is equipped with 380 GB DDR4 off-chip memory. The same Intel ICC v18 compiler was used on both platforms. We used the following compile flags: -O3, avx 512-bit flag, and -qopt-zmm-usage=high on skl. We launched the programs in native mode using the numactl -p 1 command on knl. The KMP affinity environment was set to granularity=thread,scatter for ExaStamp amr and to granularity=core,scatter for miniMD and Lammps.

Lammps is one of the most relevant classical MD simulator designed for parallel computers. The numerous built-in features for the simulation of atomic, polymeric, biological, metallic, or mesoscale systems, using a variety of force fields and boundary conditions. Lammps benchmarks include simulations up to several billion atoms on tens of thousands cores. Although the standard Lammps version is parallelized using MPI, many packages have been developed to integrate thread parallelization such as USER-OMP or USER-INTEL. The KOKKOS package has been integrated to optimize Lammps on the accelerators (Xeon Phi and GPU). A load balancing has been implemented using the shift method [START_REF] Perl | Max-min tree partitioning[END_REF] or the recursive coordinate bisection [START_REF] Marsha | A partitioning strategy for nonuniform problems on multiprocessors[END_REF] (rcb) method. Despite of specific optimizations for the knl architecture, the OpenMP version of Lammps is 10 to 15% slower than the full-MPI version [START_REF] Reinders | Intel Xeon Phi Porcessor High Performance Programing[END_REF]. Note that Lammps is compiled using the flags recommended by the documentation on Intel Xeon Phi (knl) and Xeon (skl). We have notably used the following Kokkos-runtime flags: -k on t nbThread -sf kk -pk kokkos newton on neigh half comm no. Newton's third law (neigh half) is used by all codes to take advantage of the symmetry of lj potential.

miniMD is a simplified version of Lammps written in C++. Due to its simplicity (near 5000 lines), it is possible to explore potential optimisations much more rapidly. Despite supporting only the Lennard-Jones (lj) and Embedded-Atom model (Eam) potentials, miniMD's performance and scaling behavior is representative of much larger and more complex codes with hybrid parallelization MPI/OpenMP. Simulation parameters are defined with a simple input file. miniMD is compiled using the flags recommended in the Intel Makefile, and the following runtime flags: -half_neigh=1 -t nbThread.

We tried to faithfully use the same experimental setup for min-iMD, Lammps and ExaStamp amr. Most notably, each experiment features a warm-up phase of 10 time steps, followed by 90 iterations. An amr refinement/coarsening process is required for the initialisation step. Two test cases were simulated in our experiments. The first one is a bulk material corresponding to an ideal crystal (4,000,000 atoms) with an homogeneous density. In constrat, the second one simulates a micro-jetting (3,535,584 atoms) occurring after a shock on a material. The number of atoms in these two test cases represents the quantity assigned to a node by the ExaStamp code during simulations with several hundreds of millions of atoms. The lj potential is used because its low computational cost better highlights performance issues related to suboptimal memory accesses or bus contentions (e.g. memory-bound code regions and threads synchronizations). This makes it easier to accurately identify performance issues and evaluate the effectiveness of code optimizations.

Bulk Material: a Steady-State Simulation

The first benchmark consists in simulating a bulk material featuring 4 millions of atoms in thermodynamic equilibrium. This simulation has many properties such as an homogeneous density, which also remains constant throughout the duration of the simulation. The bulk of copper is obtained by replicating an atom pattern named lattice. A Face Cubic Center (FCC) lattice of 3.54Å is used with a suitable lj potential. For each atom, 39 interactions with neighboring atoms are processed due to the 5.68Å cut-off and symmetric property of the lj potential. A radius of Verlet of 0.3Å is added to avoid re-build Verlet lists. Throughout this simulation, the tree depth D max is set to 2 to obtain best performance. A classical way to evaluate performance of MD simulations is to measure the number of atoms processed per second. The results in Table 1 report performance obtained with Lammps, miniMD and ExaStamp amr on Intel knl and skl architectures. As expected, we observe that miniMD and Lammps exhibit similar performance. ExaStamp amr is 3 times faster than miniMD, while ExaStamp amr Vec is 5.5 times faster than miniMD on the knl architecture. The difference is reduced on the skl from 3 to 1.21 and from 5.5 to 1.6 respectively. This significant difference between ExaStamp amr and the two other codes comes from the performance penalty incurred by OpenMP atomic operations used in the miniMD and Lammps kernels when updating neighboring atoms. In ExaStamp amr, the use of a wave-based task scheduling, presented in Section 3.3, avoids using such costly synchronization instructions. Moreover, the use of atomic operations in miniMD and Lammps prevents the code from being properly vectorized by the compiler. Note that performance reported for this bulk configuration, which features 4 millions of atoms, only evolves marginally when increasing the number of atoms. For instance, with 216 millions of atoms and a memory usage of about 100GB, the ExaStamp amr Vec throughput is 60 millions of atoms processed per second on the skl. The measurements show that the slight drop in performance is due to an increase of NUMA effects.

Table 2 reveals that the cost of rebuilding the amr data structures is low: under 2% on knl and 3% on skl. It also shows that the percentage of time spent in the potential computation is the lowest for ExaStamp amr, because this part has been significantly optimized: the optimization with memory blocks for vectorization reduces the time ratio of this kernel by 3-4%. Regarding the benefits gained using our cache-friendly octree traversal, we used the Linux perf profiler to finely tune the octree depth (D max ). We observed on knl that the number of cache misses of the L1 cache fluctuated according to the maximum depth of the octree. After 100 iterations, for D max =1, we observed 8.8 billions of L1 cache misses. For D max =2, the number of misses dropped to 2.6 billions, and for D max =3, it topped 9.7 billions.

On processors such as knl and skl, special attention must be paid to the vectorization of the code, because its impact on overall performance is expected to be high. Figure 8 shows the evolution of ExaStamp amr performance according to the auto-vectorization flag used at compile time (-mavx, -xMIC-AVX512, -xCore-AVX512).

The estimated speedup on the potential kernel reported by vec-report=5 option are: 2.11 with avx 512-bit knl, 1.54 with avx 512-bit skl and 1.3 with avx2 256-bit. Measurements on the full application show an acceleration of 1.88 with avx-512 on knl while it is 2.11 for the potential kernel. Similarly, we measure 1.32 with avx 512bit skl instead of 1.54 for the single kernel. The difference in the impact of vectorization on both machines is explained on the one hand by the time ratio spent in the potential kernel that is lower on the skl (see Table 2) and on the other hand because the avx 512-bit instructions set is not exactly the same.

Figure 9 describes the evolution of the simulation time of the three codes as function of the number of cores used on a processor. On both processors, the different versions of ExaStamp amr are faster than Lammps and miniMD when the number of cores is greater than 1. Note that the sequential versions of the algorithms used by Lammps and miniMD are different from the parallel versions, which explains the performance between using one core or two. In the parallel versions, there is notably the addition of OpenMP atomic instructions in the potential kernel. We can also notice in Figure 9 that on the knl, versions of ExaStamp amr have the same scalability (almost ideal) however on the skl they suffer more of NUMA effects between 24 and 48 cores than miniMD and Lammps. These effects could be reduced by using one MPI process per socket or by forcing a distribution of tasks on the sockets according to the octree position. These assumptions will constitute one of the purpose of our future developments.

Finally, this example highlights the advantages to extend ExaStamp with amr capabilities using cache blocking, task dependency graph and vectorization even for homogeneous systems. In all cases, it is faster than miniMD and Lammps. The second benchmark consists in simulating a micro-jetting [START_REF] Durand | Comparative simulations of microjetting using atomistic and continuous approaches in the presence of viscosity and surface tension[END_REF] phenomenon, which is a highly dynamic case. The micro-jetting results of a shock on a sinusoidal default in a tin material. This case performs over 3.5 million atoms with a suitable lj potential. Our experiment is calibrated so that the shock happens from the very first time steps. Thus, the density of atoms is strongly heterogeneous and the velocity of atoms is quite high. For such a highly dynamic configuration, it is expected that our amr solution behaves better than classical MD methods. Throughout this simulation, the best performance was achieved with a maximum tree depth D max equal to 2. Unfortunately, we could only compare our solution against Lammps, because miniMD is not able to cope with complex input files.

Micro-Jetting: a Highly-Dynamic Simulation

Software

Potential Table 3 reports the measurements of atoms processed per second and provides some insights into performance obtained on knl and skl processors. First of all, with an equivalent potential kernel, ExaStamp amr is 2.8 faster on knl and 1.5 faster on skl compared to Lammps. In addition, with the Vec optimization, these differences increase respectively from 2.8 to 4.7 and 1.5 to 2. As described in the previous benchmark, these differences are partly explained by a better efficiency of ExaStamp amr, especially in terms of vectorization. Throughout the simulation, 90% of cells are empty or almost empty, which have a negative performance impact on all versions except ExaStamp amr. Note that the previous generation of the ExaStamp software, which is on a par with Lammps on the bulk benchmark, suffers from a performance drop to only 1.4 million atoms per second. We have also tested the MPI version of Lammps with the rcb dynamic load balancing option. This MPI version of the code is more efficient than the OpenMP version with 6.2 million atoms per second instead of 3.3 on knl and 23.1 instead of 14.8 on skl. However, ExaStamp amr remains better in terms of CPU and memory performance (more details in Subsection 4.4). Moreover, using only one MPI process per processor using OpenMP inner parallelism will help to reduce communication costs when moving to large clusters.

Throughout this simulation, the Verlet lists are regularly updated because the atoms move in the direction of the shock. In our case, these are updated approximately every 6 time steps. As shown in the Table 4, the rebuilding of the Verlet lists costs around 20% of the simulation time for ExaStamp amr and approximately 12.5% for Lammps because it spends more time in the potential kernel and is twice slower than ExaStamp amr Vec. The time related to refinement is more significant because, at each update of the Verlet lists, the arrays of atoms within the octrees are sorted and the pointers of the cells are adjusted on these arrays. This time ratio remains proportionately cheap compared to the benefit that it provides. In addition, adapting the refinement during simulation, i. e. coarsening/refining the cells, takes around 7.5% of an iteration with construction of neighbor lists. For our case of micro-jetting, it is interesting to adapt the grid approximately every 400 time steps.

As in the bulk benchmark, we observe the impact of various compile time vectorization flags. The Figure 10 describes time with the same compilation options as used for the Figure 8. Both codes being identical, except for the potential kernel, we get the same pattern as for the Figure 8. Regarding Lammps, vectorization is more impactful thanks to the use of neighbor lists. However, this part of the code represents only 12.5% of total computation time, so it has only a moderate impact on overall performance.

A more in-depth study of the behavior of the codes is conducted in two steps, first by analyzing the times obtained at variable number of cores and second by observing the trace of the tasks with 48 cores on skl. As observed in the Figure 11, for both architectures, there is a significant difference between the sequential version and the parallel version of Lammps. As observed in the homogeneous case, ExaStamp remains much better and always shows a good adaptation to the number of cores despite the load imbalance between the different octree due to the heterogeneity of the atom distribution.

To better analyze the performance of ExaStamp amr on a large number of cores, it is worthwhile to study the task scheduling behavior induced by the wave method. Figure 12 shows a tasks execution trace on the 48 cores skl3 , using one thread per core. The trace was obtained using the VITE software [START_REF] Coulomb | Visual trace explorer (vite)[END_REF]. One can observe that only few idle slots are observed. They account for less than 5% of the simulation time. The most relevant thing to observe is that waves of tasks do not start at the same time: they overlap thanks to the use of fine-grain task dependencies. Compared to a version using a global synchronization barrier between waves, we have measured a gain of 5% on skl (48 cores) and 6% on knl (64 cores). In our implementation, thread 0 is the master thread in charge of task distribution and this is the reason why it is more often idle. At the end of the simulation, there are many threads waiting for, because there are no more ready tasks for the current iteration. This situation could be further improved by taskifying the other parts of the code, to allow the force computation step to overlap with the amr refinement, for instance.

To conclude this evaluation part, we showed that the use of octrees coupled to the wave method significantly reduces the simulation time compared to Lammps on modern processors. This is In this section we study the overall memory footprint with getrusage in Giga Bytes of both bulk and micro-jetting simulations. We expect ExaStamp amr to behave well in that respect, because one of the first expectations of amr is used to overcome the shortcomings of the linked cell method.

First, by analyzing the bulk case results in Table 5, we notice that the memory is equitably divided between the cells for ExaStamp classic and the octrees for ExaStamp amr. Regarding Lammps memory storage strategy, the atom information are stored in a structure of arrays and an array delimiting the cells with the number of atoms in each cell which minimizes memory usage. In addition to these, a storage of Verlet lists for each atom is added (more than 75% of the total storage) for each code. In the Verlet lists, miniMD and Lammps store the position of atoms while ExaStamp needs the cell index and the neighbor positions in that cell.

For the micro-jetting case, we showed the positive impact of amr on the memory footprint because a large part of the cells are mostly empty or almost empty (more than 90%). Contrary to the bulk test case, the cell storage strategy induced by the linked cell method reveals its weakness. As attested by the second column of the Table 5, octrees bear this default. Indeed, a maximum depth tree of 2 reduces by 64 the storage of empty cells. We observe that the octrees divide the memory usage by 8, even by 15 with the use of Vec option. These results are relatively close to the memory footprint of Lammps.

Thirdly, using many MPI processes, 48 in our case, on one processor can generate twice as much memory. This is due to the ghost atoms that the subdomains must communicate to each others.

amr structure drastically reduces memory footprint whenever some cells contains very few atoms, while maintaining high performance. Additionally, memory footprint may be further reduced by not storing the Verlet lists (occupying up to 75% of total memory) and recomputing them on demand. This can be achieved at the expense of a performance penalty and is required only in case of very strong memory constraints.

RELATED WORK

Many research efforts have been devoted to designing efficient MD simulations on supercomputers. Lammps and its associated miniMD mini application are probably the most widely used software in this domain, and have already been detailed previously. In the following, we carry on a short survey of significant MD software.

dl_poly_4 [START_REF] Smith | The DL POLY Classic user manual[END_REF] brings a package of subroutines, programs and data files, written in fortran 90 and designed to facilitate MD simulations and is used for production. This software is designed to performs on CPU with OpenMP and MPI, as well as on GPU.

Gromacs [START_REF] Herman | GRO-MACS: a message-passing parallel molecular dynamics implementation[END_REF] is written in C with intrinsics functions for Simd instructions and CUDA for GPU. As Lammps, Gromacs contains algorithms for neighbor search based on Verlet lists. Whereas it is designed for biochemical molecules like proteins, lipids and nucleic acids that have a lot of complicated bonded interaction, an all-atom version provides non-bonded interactions such as Lennard Jones or Buckingham potentials.

As miniMD, several simplified versions of complex MD codes have recently been developed to serve as a testbed for various optimizations over manycore architectures. The CoMD [START_REF] Mohd-Yusof | CoDesign Molecular Dynamics (CoMD) Proxy App Deep Dive[END_REF] software was designed to study the dynamical properties of various liquids and solids in the context of the ExMatEx project. The reference implementation of comd uses OpenMP, MPI, OpenCL, and can use either neighbor lists or cell-lists methods.

Moreover Lammps proposes other optimizations in the INTEL package [START_REF] Michael Brown | Optimizing legacy molecular dynamics software with directive-based offload[END_REF] such as mixed precision arithmetics, cache optimizations (tiling) and kernel offloading on knl accelerators. Optimization for GPU are also included in the Kokkos package.

CONCLUSION

Performing Molecular Dynamics N-body simulations over supercomputers is a major step towards a better understanding of physics in many scientific fields.

In this work we propose a combination of Adaptive Mesh Refinement, cache-blocking and task-based execution techniques to significantly increase the efficiency of N-body simulation on modern multicore machines, while decreasing the memory footprint. We have implemented these techniques into the ExaStamp production code, which runs on CEA's supercomputing facilities. We have evaluated the relevance of our approach by comparing performance obtained by our implementation against the state-of-the-art Lammps software over two different multicore hardware. Our experiments show that amr achieves significant performance gains over existing approaches, even for a lj potential that has low computational requirements. In the case of a dynamic, imbalanced particles configuration, ExaStamp amr can be up to 4.7 faster than Lammps on a knl processor and 2 on a skl processor. Moreover, we observed that our approach is also very effective on steady-state configurations, thanks to our cache blocking and vectorization optimizations.

Future work include the development of a hybrid MPI/OpenMP version of ExaStamp amr, focussing on a dynamic load balancing system able to move octrees between MPI processes on large-scale clusters.

Figure 1 :

 1 Figure 1: Standard methods to build atom neighbors list.

Figure 2 :

 2 Figure 2: MD simulation with amr architecture and two levels of refinement. A node is divided into octrees and the octrees are refined in leaf or inter cells.

Figure 3 :

 3 Figure3: The domain simulation (3D representation of Figure2) is decomposed into sub-domains that are distributed across MPI nodes, each sub-domain being a grid of octrees.

Figure 4 :

 4 Figure 4: Data management in a octree. The array containing all atoms information (orange color) is allocated into the Root (red) and atoms (At) are sorted by their morton index. The refined cells Inter (gray) and Leaf (blue) point on this atoms array. This strategy avoids the copy of data during refinement/coarsening steps.

  Figure 5aFigure5b

  Figure 5c

Figure 5 :

 5 Figure 5: Atom I interacts with some of the atoms inside cell (a). The Verlet lists are filled (b) by the first indexes of blocks of a vector register size (4 doubles in 256-bit) and containing a minimum of one atom in the radius of Verlet. During the potential step (c), the atoms are loaded by blocks and could be pre-fetched.

Figure 6 :

 6 Figure 6: 2D representation of a part of the MD simulation. The atoms are distributed into the cells themselves included into octrees with a tree depth of 1. The octrees are colored in 4 colors such that two atoms included in two octrees of the same color do not interact together.

Figure 7 :

 7 Figure 7: Task dependency graph in 2D. octree are classified into 4 waves (8 in 3D), cells classified in the wave (n) n >=0 are unlocked by specific cells in the wave (n -1) n >=1 .

Figure 8 :

 8 Figure 8: Time simulation for different vectorization flags on Ex-aStamp amr Vec, the version no Vec corresponding to ExaStamp with avx-512 bits and no option Vec. Software performed on 64 cores on knl (left) and 48 cores on skl (right).

Figure 9 :

 9 Figure 9: Simulation time (lower is better) with log axis for miniMD, Lammps, ExaStamp amr and ExaStamp amr Vec. For both architectures, the codes keep the same ranking regarding simulation time, with ExaStamp amr Vec being the fastest, then ExaStamp amr, then miniMD and finally Lammps.

Figure 10 :

 10 Figure 10: Simulation time for different vectorization flags on Ex-aStamp amr Vec, the version no Vec corresponding to ExaStamp with avx-512 bits and no option Vec. Software performed on 64 cores on knl (left) and 48 cores on skl (right).

Figure 12 :

 12 Figure 12: Visualization of tasks execution throughout a time step of a micro-jetting simulation on skl on 48 cores. The colors represent the type of task, from light green to light blue through black corresponding to the tasks from wave 1 to 8. The red portions correspond to the moment when a thread is idle.

Table 2 :

 2 Time ratio measured for different parts of the bulk of copper simulation: potential kernel (Potential) and related to the amr structure (Refinement). Time ratio of other parts are not reported (Scheme, Reorganisation, Update Ghost). Note that the Verlet lists are not updated.

	Software	Potential	Refinement
	Lammps knl	96.4%	NA
	Lammps skl	88.3%	NA
	miniMD knl	97.5%	NA
	miniMD skl	91.2%	NA
	ExaStamp amr knl	94.8%	0.9%
	ExaStamp amr skl	81.1%	1.9%
	ExaStamp amr Vec knl	91%	1.7%
	ExaStamp amr Vec skl	77.3%	2.6%

Table 3 :

 3 Number of atoms processed per second (higher is better) by Lammps, ExaStamp amr with and without Vec option (mean).

	Simulation code	million atoms processed per second
	Architecture	Knights Landing	Skylake
	Lammps	3.3	14.8
	ExaStamp amr	9.2	22.3
	ExaStamp amr Vec	15.5	29.1
	Software performed on 64 cores on knl (left) and 48 cores on skl
	(right).

Table 4 :

 4 Time ratio measured for different parts of the microjetting simulation: potential kernel (Potential), building Verlet lists (Verlet Lists) and related to the amr structure (Refinement). Time ratio of other parts are not reported (Scheme, Reorganisation, Update Ghost.)

	Verlet	Refine-
	Lists	ment

  the simulation

					ExaStamp amr KNL					ExaStamp amr SKL
		10 3			ExaStamp amr Vec KNL MiniMD KNL				ExaStamp amr Vec SKL MiniMD SKL
					Lammps kokkos KNL	10 2				Lammps kokkos SKL
	Time(s)	10 2					Time(s)
								10 1
		10 1
		1	2	8	16	32	64	1	2	6	12	24	48
				Cores						Cores

Table 5 :

 5 Figure 11: Simulation time (lower is better) of a micro-jetting with log axis for Lammps, ExaStamp amr and ExaStamp amr Vec. For both architectures, the codes keep the same ranking regarding simulation time, with ExaStamp amr Vec being the fastest, then ExaStamp amr, then miniMD and finally Lammps. Memory footprint (lower is better) of the different codes for both datasets: Bulk and Micro-jetting.

				ExaStamp amr KNL					ExaStamp amr SKL
				ExaStamp amr Vec KNL				ExaStamp amr Vec SKL
	10 3			Lammps kokkos KNL					Lammps kokkos SKL
	Time(s)					Time(s)	10 2
	10 2
							10 1
	1	2	8	16	32	64	1	2	6	12	24	48
			Cores						Cores
	explained by the good use of cache blocking techniques, vectoriza-
	tion and task based scheduling.
	4.4 Memory footprint
	Simulation code		Memory foot print (GB)
	Dataset		Bulk	Micro-jetting
	Lammps		2.3		6.4
	Lammps (48 MPI processes)	3.4		13.8
	miniMD		2.1		NA
	ExaStamp Classic		3.3		48
	ExaStamp amr		3.7		8
	ExaStamp amr Vec		1.9		3.2

The analysis presented focuses on the use of the skl with

cores but the comments are also relevant for the knl.

ACKNOWLEDGMENT

This work was funded by the French Programme d'Investissements d'Avenir (PIA) project SMICE.


Tristan Monchal

Amina Ndiaye

Blandine Gadegbeku

Etienne Javouhey

Olivier Monneuse

Abdominopelvic injuries due to road traffic accidents: Characteristics in a registry of , victims

Keywords: Road traffic accident, abdominal trauma, pelvic trauma, trauma registry, Abbreviated Injury Scale

Objective: Road traffic accidents (RTAs) are the first cause of abdominopelvic injuries (APIs). The objective of this study was to describe the characteristics and severity of APIs due to traffic accidents in a large French trauma registry and to identify risk factors for API. Methods: All victims from the French Rhône registry of victims of RTAs were analyzed from 1996 to 2013. This registry contained data that were issued over a 20-year period from 245 medical departments, from prehospital care until re-adaptation, and forensic medicine departments. All APIs, defined as an injury between the diaphragm and the pelvic bone, were extracted and studied. Results: Among 162,695 victims, 10,165 had an API (6.7%). Accidents frequently involved young men and 2 cars. Mean Injury Severity Score (ISS) was 8.7. Mortality rate was 5.6%. Soft tissue injuries largely predominated (n = 6,388; 54.4% of patients). Overall, 2,322 victims had a pelvic bone injury. Internal abdominal organs were involved in 2,425 patients; the most frequent were the spleen, liver, and kidney. Wearing of the seat belt appeared to be a significant protective factor in API, including serious injuries. A partial analysis over the past 2 years among the most severe patients hospitalized in the intensive care unit indicated that nonoperative management was carried out in two thirds of the wounded. In uni-or multivariate analysis, sex, age, type of user, antagonist, time of occurrence, associated severe lesions, or wearing of the seat belt were statistically associated with the occurrence of API, highlighting a more dangerous user profile. Conclusions: Abdominopelvic injuries concern a minority of road traffic injuries, but they are responsible for significant mortality. Large solid organs are the most frequently affected. Women drivers wearing a seat belt and driving in town during the day appear to be more protected against API.

Introduction

Road traffic accidents (RTAs) are still the first cause of abdominopelvic injuries (APIs; [START_REF] Moore | Abdominal and pelvic injuries[END_REF][START_REF] Smith | Abdominal trauma: a disease in evolution[END_REF]. They are responsible for important mortality, which makes them a public health problem [START_REF] Ndiaye | The fatal injuries of car drivers[END_REF][START_REF] Tournier | A few seconds to have an accident, a long time to recover: consequences for road accident victims from the ESPARR cohort 2 years after the accident[END_REF]. Few series in the literature have attempted to describe the epidemiology of APIs, particularly related to RTAs [START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF][START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF][START_REF] Helmi | Abdominal trauma due to road traffic accidents in Qatar[END_REF].

For more than 20 years, the Rhone Registry of Victims from Road Traffic Accidents (Registre des Victimes d' Accidents de la Route du Rhône = ReVARRhône) has been collecting exhaustive data from all medical facilities in the French Rhone department on victims of RTAs [START_REF] Gadegbeku | Victimes d'accidents de la circulation: quelles évolutions? Registre du Rhône, France, 1996-2008[END_REF].

The objective of this study was to describe the frequency and characteristics of abdominopelvic trauma from this registry and to study the risk factors among several available variables.

CONTACT Tristan Monchal tristan.monchal@free.fr Department of General Surgery, Sainte Anne Military Hospital, BP ,  Toulon cedex , France. Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/gcpi. Associate Editor Sjaan Koppel oversaw the review of this article.

Methods

Overview of the Rhône Registry

The Rhone Registry, certified by the French National Committee of Registries, is a unique model in Europe; since 1995 it has been collecting data pertaining to the victims of road traffic accidents [START_REF] Gadegbeku | Victimes d'accidents de la circulation: quelles évolutions? Registre du Rhône, France, 1996-2008[END_REF]. It is part of the French Institute of Transport, Development and Road Network Science and Technology (Institut Français des Sciences et Technologies des Transports, de l' Aménagement et des Réseaux, IFSTTAR). Its function is to perform the continuous recording of all physically injured victims due to RTAs involving at least one moving vehicle (including skates and skateboards) and occurring in the geographic area of the Rhône department, with a population of 1,779,845 inhabitants (2013 data). The majority of the victims were residents of the department.

The collection is based on 245 public and private health facilities dealing with trauma care: Fire and rescue aid services, prehospital care services, emergency departments, displacement, resuscitation and intensive care, surgical and medical departments, forensic medicine, and rehabilitation, including in neighboring departments. Various types of information are collected regarding the patient, the accident, the lesions, and outcomes. Each victim has a unique identification number. Victims or their families are asked whether it is necessary to supplement some missing data. The registry team is in charge of centralization, verification of the data from different sources about the same accident or victim, coding, storage and filing, and statistical analysis. The Abbreviated Injury Scale (AIS 90) was used to code all injuries (Association for the Advancement of Automotive Medicine 1990). The Maximum Abbreviated Injury Scale (MAIS) is the highest (maximum) immediate gravity score among all the injuries. The aims of this registry, used by the public services, are to improve knowledge of traffic accidents and to help decrease the number and gravity of these accidents.

Study design

We conducted an observational study from January 1 1996 to December 31, 2013. During this time frame, all victims with an abdominopelvic injury, including diaphragm and pelvic bone, were analyzed (AIS codes R = 5 for the abdomen, 4.4.06 for the diaphragm, and 8.5.26 + 8.5.28 + 8.5.30 for the pelvic bone). Several data were extracted and analyzed: Demographics (sex, age), characteristics of the accident (type of user, protection devices, antagonist, type of road, time), anatomical lesions, and elements of gravity (Injury Severity Score [ISS], AIS, mortality).

Statistical analysis

Statistical analyses were performed using Statistical Analysis System 9.4 (SAS Institute Inc., Cary, NC). A P value of .05 was used to define statistical significance. Chi-square test, or Fisher's exact test when the cell sizes were not large enough, was used to compare the equality of proportions. Logistic regression analysis was used to explain the presence of an API depending on different factors.

Results

Overall characteristics

During the time frame of the study, there were 204,227 forms collected, which were related to 143,027 accidents and a total of 162,695 victims. Of these victims, 10,165 presented with APIs (6.2%). Mean annual incidence of APIs was 565 victims (34/100,000 inhabitants in Rhône).

The mean age was 32.2 years; those aged 20-24 years were most frequently involved. Overall, 62.4% of APIs involved men (6,347 men and 3,815 women, ratio 1.7).

Almost half of the victims were car drivers (46%), and around one quarter were motorcyclists (24.7%), followed by pedestrians (12.8%), bicyclists (11%), and rarely other vehicles. The most common antagonist was a car (55.4%), followed by no antagonist (22.4%), fixed obstacle (9.3%), and truck or bus (8.2%).

Severe APIs (AIS 3+) involved 2,676 victims (26.3% of APIs and 1.6% of the whole registry). The mean ISS was 8.7 (24.6 for the AIS 3+ subgroup). The mortality rate for APIs was 5.6% (1.1% in the entire registry), with a majority of deaths occurring on scene (314 patients, 55.3%) or during the first 24 h (118 patients, 20.8%).

Concomitant injured body areas, frequencies, and AIS gravity are presented in Table 1. The mean number of injuries was 2.5 per victim. The most frequently extra-abdominal-associated injured zones included lower limbs (with the exception of the pelvic bone), chest (with the exception of the diaphragm), upper limbs, head, and spine. Among the most severe AIS 4+ patients, chest was the most frequently involved (41%), followed by head (27%) and abdomen (21%). We found the same frequency for fatalities in our study when data were available due to a forensic examination or a postmortem computed tomography (CT) scan. In patients with both severe AIS 3+ cerebral and AIS 3+ abdominopelvic injuries, the mortality rate increased to 65%.

The majority of lesions were superficial ones, involving skin and soft tissues (n = 6,388; 54.4%). There were 2,322 pelvic bone injuries and 72 diaphragm injuries. Among the intra-abdominal lesions (n = 2,497; 21.3%), the large solid organs were the most commonly injured: spleen, liver, and kidneys (Figure 1). Road user type had no significant impact on the distribution of most of the injured organs. On the other hand, injuries to external genitalia (n = 462; 3.9%) occurred almost exclusively among 2-wheeler drivers.

Impact of wearing a seat belt among car drivers

Among the car users (n = 4,675), 75.1% wore a seat belt, 12.1% did not, and this information was missing in 12.8%. Wearing a seat belt appeared to be protective against severe lesions: AIS 2 to 5 injuries and mortality were significantly higher among nonbelted victims.

Distribution of the injuries showed some differences depending on whether or not a seat belt was worn: Splenic, renal, pancreatic, adrenal, diaphragmatic, and retroperitoneal injuries were more frequent among nonbelted users, whereas differences for the other organs, including liver or hollow viscous injuries, were not significant.

Hospitalization in the intensive care unit for severe APIs

Among AIS 3+ severe APIs (n = 2,676), 1,537 people (57.4%) were admitted to the intensive care unit. During the past 2 years of our study (2012)(2013), analysis of medical reports from the intensive care unit for AIS 3+ injuries (n = 181) revealed that 38 patients (21%) required a hemostatic surgical procedure, 19 (10.5%) had a radiological arterial embolization, and the majority (124, 68.5%) could be treated with nonoperative management.

Analysis of risk factors for APIs

A multivariate model was used to ensure that the interaction between all factors was taken into account (for example, motorcyclists are mostly men, whereas women are more often belted [START_REF] Gadegbeku | Séquelles majeures en traumatologie routiére, registre du Rhône, 1996-2003 [Major incapacity after roadtraffic traumatisms, Rhône registry, France, 1996-2003[END_REF]; Table 2). Victims of APIs were more often men (odds ratio [OR] = 1.111; P = .005). Compared to the reference group of 26-to 35-year-olds, children <15 years old were less affected (OR = 0.835; P = .006); whereas those older than 46 years old were increasingly more affected (P < .001 for each). Pedestrians and 2-wheeler users had greater risk compared to car drivers (OR > 2 with P < .001). The absence of obstacles was a protective factor (OR = 0.638; P < .001), whereas a fixed obstacle or a heavy goods vehicle was at higher risk than a car for antagonists (P < .001). Regarding the road type, compared to a city street, the risk of having an API was higher on a road (OR = 2.174; P < .001), not significant for a highway (P = .253), and lower for all other types such as a path, a parking area, etc. (OR = 0.838; P < .001). Analysis of the time of the accident revealed that there was a lower risk of API between 8 a.m. and 3 p.m.; on the other hand, the negative impact of nighttime from 8 p.m. to 7 a.m. was clear (P < .001).

A severe AIS 3+ injury to the thorax (excluding diaphragm), head (AIS 3+), or spine (AIS 2+) was a major risk factor of associated with API (all OR > 10; P < .001). Finally, not wearing a seat belt was also a significant risk factor (OR = 2.556; P < .001).

Discussion

Our study reveals a 6.2% incidence of abdominopelvic injuries within the whole registry. Other series in the literature reported higher rates, between 11 and 18%, but were either limited to inhospital series or included other causes of trauma, making them difficult to compare (American College of Surgeons Committee on Trauma [ACSCOT] 2016; [START_REF] Arumugam | Frequency, causes and pattern of abdominal trauma: a 4-year descriptive analysis[END_REF][START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF][START_REF] Helmi | Abdominal trauma due to road traffic accidents in Qatar[END_REF]. [START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF] reported, in a series that also included prehospital care interventions, a 9.3% incidence. The Rhône Registry also collects data pertaining to untransported self-consulting patients, frequently with minor lesions. Moreover, our study commenced 20 years ago, radiological exams were not well organized, and abdominal lesions were probably often underdiagnosed, whereas the current systematization of total body CT scan allows a better detection and characterization of injuries. Finally, the Rhône department is a small and urban area, where the gravity of RTAs is less severe than in the rest of the country.

The incidence of APIs in the recent years of our study decreased to 24.5/100,000 inhabitants in 2013. At the same time, the overall incidence of RTAs has been continuously decreasing, due to French laws regarding speeding and driving under the influence of alcohol from the early 2000s, combined with advances in car engineering and road safety, similar to other countries [START_REF] Ernstberger | Decrease of morbidity in road traffic accidents in a high income country-an analysis of 24,405 accidents in a 21 year period[END_REF].

The predominance of young men constitutes constant data in road traffic accident registries, even if a shift toward an increase in the age of victims can be observed (ACSCOT 2016;[START_REF] Lapostolle | The burden of road traffic accidents in a French departement: the description of the injuries and recent changes[END_REF].

The ISS was used to describe the severity of injuries. The new ISS (NISS), using the 3 most severe injuries regardless of body region, is a better predictor of survival and could have been adapted here to take into account the association of multiple abdominal injuries [START_REF] Brenneman | Measuring injury severity: time for a change?[END_REF]. According to the period of data collection, the ISS was the most easily available scale.

Abdominopelvic injuries are quite rare and also rarely severe, but our study confirms that they remain associated with a very significant mortality: 5.6% vs. 1.1% in the rest of the registry (P < .0001; [START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF]. Mortality occurs mostly on scene and during the "golden hour, " which can be explained by uncontrollable hemorrhaging and severe head trauma. Only 10% of the mortality is related to APIs, after thoracic and head injuries. The association of severe cerebral and abdominopelvic injuries was associated with a mortality rate of 65%. This confirms the well-known deleterious combination of severe hemorrhage, such as severe abdominal injury, and a cranioencephalic injury [START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF][START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF]. However, in this study, it was not possible to systematically identify the reason for death, which was mainly due either to cerebral injuries or to massive bleeding from the chest, abdomen, or pelvis. Because road traffic accidents are considered to be a cause of death in their own right, autopsies are rarely performed in France for such fatalities [START_REF] Ndiaye | The fatal injuries of car drivers[END_REF]). In those cases, anatomical injuries could be identified with an external physical examination only or a virtual forensic exam with a CT scan in hospitals where this protocol existed. Soft tissue injuries, poorly studied because of their frequent benignity, account for more than half of lesions. The pelvic bone, as reported by [START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF], is the most frequently injured organ among APIs. Its overall incidence within the whole registry is 1.4%. Due to high-velocity trauma, which poses a high risk of associated lesions, fractures, or disruptions of the pelvic bone, is responsible for increased mortality [START_REF] Sathy | The effect of pelvic fracture on mortality after trauma: an analysis of 63,000 trauma patients[END_REF]. As expected for blunt abdominal trauma, solid organs are the most frequently injured intra-abdominal organs: First the spleen and liver in similar proportions, followed by the kidneys [START_REF] Arumugam | Frequency, causes and pattern of abdominal trauma: a 4-year descriptive analysis[END_REF][START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF][START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF][START_REF] Helmi | Abdominal trauma due to road traffic accidents in Qatar[END_REF]. External genitalia are relatively frequently involved (3.9% of all APIs), mostly among 2-wheeler drivers, and remain anecdotic among other drivers.

The seat belt is a protection device that contributes to decreasing the severity of road traffic accidents, and this study logically confirms its beneficial effect [START_REF] Gadegbeku | Victimes d'accidents de la circulation: quelles évolutions? Registre du Rhône, France, 1996-2008[END_REF][START_REF] Nash | Seat belt use and its effect on abdominal trauma: a national trauma databank study[END_REF]. Not wearing a seat belt is not only an independent risk factor of having an API in the multivariate model, but, among victims of API, the seat belt is a protective factor against severe lesions and mortality. We observe a significant decrease in most internal solid organ injuries due to the absence of a seat belt but there is no increase in hollow viscus injuries, although commonly described in the literature [START_REF] Nash | Seat belt use and its effect on abdominal trauma: a national trauma databank study[END_REF][START_REF] Sharma | Clinical implications of the seat belt sign in blunt trauma[END_REF]. Seat belt marks are not specifically documented in our registry but are anatomically included with the other abdominal cutaneous hematomas of the AIS scale. Thus, we could not study the relation between seat belt mark as such and APIs but only the relation between being restrained (with or without a seat belt sign) and APIs. This is reported as a sign of gravity, related to the violence of the accident, but it is not specific to hollow viscus injuries (Velhamos et al. 1999).

In multivariate analysis, all factors in our analysis were significantly correlated with the risk of having an API. The type of vehicle or locomotion is the most determinant factor: 2-Wheelers and pedestrians are more exposed, probably due to less protection. No third party is a protective factor, in contrast to the kinetics added by an antagonist. A growing vulnerability with age also appears [START_REF] Brasel | What mechanism justifies abdominal evaluation in motor vehicle crashes?[END_REF]. The impact of place and time results in a maximum risk on a rural road, as well as at night, and may be related to the violence of impacts under possible conditions of high speed, frontal antagonists, and low visibility. The specific high risk linked to the male sex may be explained by a higher average speed or a greater risk behaviors (alcohol or narcotics); these factors are unfortunately not available in the registry [START_REF] Martin | Différences entre les hommes et les femmes face au risque routier [Differences between males and females in road-traffic accidents[END_REF]). According to our study, the user who would be most protected against API is a belted woman driving a car during daytime on a city street, whereas an elderly pedestrian or a motorized 2-wheeler hitting a heavy vehicle at night on a rural road appears to be much more exposed to the occurrence of an API, keeping in mind that this is only a statistical association but that these might be proxy variables for other unrecorded variables, like speed, for instance.

For more than 20 years, the Rhône Registry has provided an exhaustive census of all victims of a traumatic injury after a traffic accident involving a means of locomotion and occurring on a road open to traffic. To our knowledge, this unique model has no equivalent at a national level or in Europe. In the United States, the National Trauma Data Bank is a national registry that is administered by the American College of Surgeons Committee on Trauma (ACSCOT 2016). It is the largest cohort in the world, with more than 7 million registered casualties (more than 861,000 inclusions in 2015), with data coming from 747 U.S. and Canadian volunteer trauma centers. The Rhône Registry does not suffer from hospital admission bias because it is based on the totality of the department's emergency and care services. It aims at completeness in both public and private hospital networks [START_REF] Gadegbeku | Séquelles majeures en traumatologie routiére, registre du Rhône, 1996-2003 [Major incapacity after roadtraffic traumatisms, Rhône registry, France, 1996-2003[END_REF]. However, some victims may not be included. Completeness is estimated at 73% for minor and moderately wounded survivors and at 87% for seriously wounded survivors [START_REF] Amoros | Estimating non-fatal road casualties in a large French county, using the capture-recapture method[END_REF]. Police statistics, which are the basis of road safety studies, underestimate the number of victims injured in traffic accidents [START_REF] Amoros | Under-reporting of road crash casualties in France[END_REF]). In addition, victims who directly consult a physician or take care of themselves are not recorded in the registry.

The information provided by the registry has some limitations. First, due to the medical source, it is not possible to assess detailed data about mechanisms such as velocity, rollover, impact direction, etc. Velocity and rollover were shown to be independently associated with the presence of abdominal injury [START_REF] Brasel | What mechanism justifies abdominal evaluation in motor vehicle crashes?[END_REF], whereas the side of impact might induce significant differences in the type of injuries and outcomes [START_REF] Charyk | A comparison of injuries, crashes, and outcomes for pediatric rear occupants in traffic motor vehicle collisions[END_REF]. Alcohol, as well as other drugs, is another significant factor in RTAs that has been an important area for injury prevention efforts [START_REF] Missoni | Alcohol-related road traffic accidents before and after passing of the Road Traffic Safety Act in Croatia[END_REF]. These data are unfortunately not available in the registry but have been the subject of several parallel studies by our team [START_REF] Martin | Cannabis, alcohol and fatal road accidents[END_REF]. Lastly, the main limitation may be the lack of directly available data regarding physiology and the management of trauma, which could have been informative and potentially relevant with severe hemorrhagic APIs [START_REF] Cheynel | Abdominal and pelvic injuries caused by road traffic accidents: characteristics and outcomes in a French cohort of 2,009 casualties[END_REF][START_REF] Haddad | Profile, outcome and predictors of mortality of abdomino-pelvic trauma patients in a tertiary intensive care unit in Saudi Arabia[END_REF]. It was still possible to perform a partial analysis over the past 2 years, by searching records of patients with severe AIS 3+ APIs and admitted to an intensive care unit. This confirmed the importance of nonoperative management even for severe lesions, due to progress in radiology and intensive care.

However, the anatomical descriptions are based on information from several departments, enabling very precise and fine data aggregation. The size, completeness, and accuracy of the registry make it a reliable epidemiological tool, which is used by governments or industries to improve knowledge and the prevention of road traffic accidents.

Our study confirms that the incidence of APIs in RTAs is low, but that they are responsible for an increased mortality. The pelvic bone is frequently affected in more than one fifth of cases; spleen and liver are the most frequently injured organs, each one in more than one quarter of the cases. We did not find any increased risk of hollow viscus injuries due to the seat belt, which appears to be a significant protective factor for severe lesions. Our statistical analysis reveals that a belted female car driver during daytime in town could be the best protected against abdominopelvic injuries. Further information about physiology or treatment would be relevant in the future for a better comprehensive clinical approach to those injuries.

Table  .

  Concomitant injuries with ISS according to MAIS.

	27.2%
	25.2%
	25%
	20%
		15.8%
	15%
	10%		7.8% 7.8%
			6.0%
	5%			4.4%	2.9% 2.8%
	0%
	Body region	MAIS = 	MAIS = 	MAIS = 	MAIS = 	MAIS = 	MAIS = 	Total (AIS +)
	Head								, ()
	Face	,							, ()
	Neck								 ()
	Chest	,							, ()
	Abdomen	,							, ()
	Spine	,							, ()
	Upper limb	,							, ()
	Lower limb	,	,		,				, ()
	Skin								 ()

30%

Figure . Overall distribution of internal abdominal injuries (n = , victims).

Table  .

  Risk factors for abdominopelvic injuries: Odds ratios estimated from logistic regression including , victims of traffic accidents between  and .

			Univariate analysis			Multivariate analysis
		OR	% Confidence interval	P value	AOR	% Confidence interval	P value
	Sex
	Male	.	[.-.]	<.	.	[.-.]	.
	Female		-	-		-	-
	Age (years)
	-	.	[.-.]	.	.	[.-.]	.
	-	.	[.-.]	.	.	[.-.]	.
	-		-	-		-	-
	-	.	[.-.]	.	.	[.-.]	.
	-	.	[.-.]	<.	.	[.-.]	<.
	-	.	[.-.]	<.	.	[.-.]	<.
	-	.	[.-.]	<.	.	[.-.]	<.
	>	.	[.-.]	<.	.	[.-.]	<.
	Type of user
	Pedestrian	.	[.-.]	<.	.	[.-.]	<.
	Car		-	-	-	-	-
	Motorcycle	.	[.-.]	<.	.	[.-.]	<.
	Bicycle	.	[.-.]	.	.	[.-.]	<.
	Antagonist
	None	.	[.-.]	<.	.	[.-.]	<.
	Car		-	-	-	-	-
	Truck/bus	.	[.-.]	<.	.	[.-.]	<.
	Fixed obstacle	.	[.-.]	<.	.	[.-.]	<.
	Other	.	[.-.]	.	.	[.-.]	.
	Road network
	Highway	.	[.-.]	.	.	[.-.]	.
	Rural road	.	[.-.]	<.	.	[.-.]	<.
	City street		-	-	-	-	-
	Other	.	[.-.]	<.	.	[.-.]	<.
	Time of accident
	:-: a.m.	.	[.-.]	<.	.	[.-.]	<.
	:-: a.m.	.	[.-.]	<.	.	[.-.]	<.
	:-: a.m.	.	[.-.]	.	.	[.-.]	<.
	:-: p.m.	.	[.-.]	.	.	[.-.]	.
	:-: p.m.		-	-	-	-	-
	:-: p.m.	.	[.-.]	<.	.	[.-.]	<.

Note: OR = odd ratio; AOR = adjusted odd ratio.

©  Taylor & Francis Group, LLC

Acknowledgments

The authors thank all of the people who took part in data collection and filing for the Association for the Registry of Victims of Road-Traffic Accidents in Rhône (ARVAC) and the French Institute of Science and Technology for Transport, Spatial Planning, Development and Networks (IFSTTAR-UMRESTTE).


Charles Paillet
email: paillet@lmt.ens-cachan.fr

Pierre Ladevèze

David Néron

Parameter-Multiscale PGD Methods for High Dimensional Parametric Spaces

Keywords: Reduced Order Model (ROM), PGD, Multiparametric

come   L'archive ouverte pluridisciplinaire