dataesr/scientific-paragraphs-categorization huggingface dataset:
This article describes how we manage to increase performance and to extend features of a large parallel application through the use of simultaneous multithreading (SMT) and by designing a robust parallel transpose algorithm. The semi-Lagrangian code Gysela typically performs large physics simulations using a few thousands of cores, between 1k cores up to 16k on x86-based clusters. However, simulations with finer resolutions and with kinetic electrons increase those needs by a huge factor, providing a good example of applications requiring Exascale machines. To improve Gysela compute times, we take advantage of efficient SMT implementations available on recent INTEL architectures. We also analyze the cost of a transposition communication scheme that involves a large number of cores in our case.
INTRODUCTION A key factor that determines the performance of magnetic plasma containment devices as potential fusion reactors is the transport of heat, particles, and momentum. To study turbulent transport and to model Tokamak fusion plasmas, several parallel codes have been designed over the years. Computational resources available nowadays have allowed the development of several petascale codes based on the wellestablished gyrokinetic framework. In this article, we focus on the gysela gyrokinetic code parallelized using a hybrid MPI/OpenMP paradigm  [1, 14, 15, 18]. The main data in gysela is an ion distribution function in a 5D phase space. (r, θ, ϕ) are the space variables with r the radial direction, θ and ϕ respectively the poloidal and toroidal directions and (v , µ) the velocity variables.
v is the velocity parallel to the magnetic field lines while the magnetic momentum µ is proportional to the perpendicular velocity squared. The parallelization of the code take advantage of the fact that it is an adiabatic invariant therefore plays the role of a parameter. Uniform mesh partition is used, the size of the phase space grid is Nr N θ Nϕ Nv Nµ. The computational domain is defined on (r, θ, ϕ, v , µ) ∈ [rmin, rmax] × [0, 2π] × [0, 2π] × [vmin, vmax] × [µmin, µmax]. Lower and upper bounds in r, v and µ are fixed according to the physical case. One of the main computational cost is the advection operator (that performs a single time step integration of the Vlasov equation). This operator is split in order to avoid a costly 4D advection in (r, θ, ϕ, v ). The chosen Strang splitting consists in four 1D advections (along ϕ and v directions) and one 2D advection (both r and θ directions are treated simultaneously), that are applied at each time step  [7, 8].
GYSELA APPLICATION The gysela code is a non-linear 5D global gyrokinetic full-f code which performs flux-driven simulations of ion temperature gradient driven turbulence (ITG) in the electrostatic limit with adiabatic electrons. It solves the standard gyrokinetic equation for the full-f distribution function, i. e. no assumption on scale separation between equilibrium and perturbations is done. This 5D equation is self-consistently coupled to a 3D quasineutrality equation. The code also includes other features not described here (ion-ion collisions, several kind of heat sources). The code has the originality to be based on a semi-Lagrangian scheme  [20]  and it is parallelized using an hybrid OpenMP/MPI paradigm  [5, 14]. The presence of a strong magnetic field in the Tokamak induces a large anisotropy of the transport along and across the magnetic field lines.
Therefore, switching from Algo. 1 to Algo. 2 requires to adapt the communication scheme to send the appropriate J0 φderivatives to each MPI process. In Algo 1, the domain decomposition of these derivatives matches the main code decomposition (r = local, θ = local, ϕ = * , v = * , µ = local), they are used equally into 1D and 2D advections. However, in Algo. 2, the new transposed 2D advection algorithm needs additional communications to get the J0 φ-derivatives on the new domain decomposition (r = * , θ = * , ϕ = local, v = local, µ = local). The costs of these additional communications will be exhibited in the next subsection.
Experimental Evaluation Timing measurements of this subsection have been realized on the Curie machine at the French GENCI-TGCC-CEA computing center. Each computing node is a dual socket Intel Xeon E5-2680 (Sandy Bridge): 2×(8 cores, 2. 7 GHz, 20 MiB shared L3-cache, DDR3 1600 Mhz memory) with 64 GiB of RAM. Let us mention that the thread level for MPI is set to MPI THREAD FUNNELED in gysela. MPI THREAD MULTIPLE approach is avoided because of portability issues on production machines (notably with BullxMPI proposed on several supercomputers). Furthermore, it is not uncommon that codes relying on the level of MPI THREAD MULTIPLE run slower than the case where one of the other modes has been chosen  [24]  (management of locks and thread-safety can be costly). For a strong scaling experiment, we choose a domain size of Nr × N θ × Nϕ × Nv × Nµ = 512 × 512 × 128 × 128 × 32, representing 1 TiB of data for a single distribution function.
The relative efficiencies of Fig. 3  and 4  for the two algorithms are very similar. For the entire application this efficiency is competitive and larger than 60% at 65536 cores in this strong scaling benchmark. Practically, let us notice that physicists run the code between 1k up to 8k cores commonly for this type of domain size (16k cores is uncommon). Typically, we avoid the cases with 32k and 64k cores presented in Fig. 3  and 4. Thus, we mainly target high parallel efficiency for production runs in order to maximize the number of results one can get within the amount of CPU hours that we obtain every year on several supercomputers. On Fig. 5 , 6 and 7, a zoom on the time measurement of the new method compared to the original algorithm is given. Both computations and communications are accounted for in the curves (rationale is that separation of computation versus communications measurements is not easy for subroutines that finely mix both steps).
The transpose-based algorithm is one example of optimization that is favorable in terms of computation at the expanse of network bandwidth use. This statement should however be balanced by two observations. First, we can cope with a tripling of this relative cost without incurring a severe penalty on the total execution time. Second, the main problem expected regarding Exascale networks is related to latency, it is not so clear for bandwidth. In any case, we will have to keep comparing the halo vs. transpose based approaches with each new hardware architecture. If the situation was to become too critical, we could consider introducing a pipelined version of the transpose communication pattern in order to partially overlap communication costs with the upcoming computations. Using the current programming model to express such a pipeline might however make the code difficult to read and switching to a task-based model where computations are automatically scheduled when the data becomes available would likely help in this case.
This enables to leverage both units even for applications not perfectly balanced in term of multiply/add. The addition of new execution units combined with the improved capabilities of the existing ones means that it becomes more realistic for the scheduler to submit close to its maximum of 4 µ-operations per cycle thus improving the parallelism at this level. As with any kind of parallelism, the difficulty with this design is for the code to expose enough in-fly µ-operations from which the scheduler can choose. This is why simultaneous multi-threading 1 (SMT) is more and more important Figure  9 : Haswell CPU Core Pipeline Functionality, extract from  [10]  to feed all units (also known as a reduction of horizontal waste). In addition, using multiple hardware threads makes it possible to hide latencies related to data access by switching to another thread that is not waiting for data (reduction of vertical waste).
Performance Metrics From the user's point of view, the impact of SMT is typically measured by calculating the relative speed-up attained, e. g. the code is sped up by x% using SMT compared to one thread per core configuration. In order to measure the core-level effects of SMT, a useful quantity to analyze is the utilization level of the core's execution units. In particular, it is of utmost importance to access a significant fraction of micro-operation slots (executing units) that are available to execute an instruction. If the number of instructions performed per cycle (IPC) is high, then the execution units are being kept busy doing useful work. We will measure IPC in some important kernels of gysela. Let us notice that on Haswell processor, IPC can reach a maximal value of 4. This is due to instruction retirement and decode units that can treat up to 4 micro-operations per cycle  [9].
As a result, the simultaneous multi-threading with 2 threads per core gives a benefit of 21% up to 28% over the standard execution time (deployment with one thread per core). While an improvement is expected with SMT, as already reported for other applications (e. g. [11, 23] ), this speedup is quite high for a HPC application. We have investigated the most intensive computation parts of the code with Paraver tools (www. bsc. es/paraver). The tools are based on traces capturing the detailed behavior of the different MPI processes and threads along time. Calls to the MPI and OpenMP runtime can be enriched with hardware counters, so we were able to measure the instructions and cycles for each computation region. We observe that for each intensive computation kernel the number of instructions per cycle (IPC) cumulated over the 2 threads on one core with SMT is always higher than the IPC obtained with one thread per core without SMT.
See for example 2D advection and Transpose steps (Useful duration plots), there is much black color at the end of these steps. This is due to final MPI barriers. Nevertheless the impact is relatively low in this reduced test case because the tool reported a parallel efficiency of 97% over the entire application indicating that only 3% of the iteration time is spend on the MPI and OpenMP parallel runtimes. The impact is stronger on larger cases, because load imbalance is bigger. 3. The transpose steps show a lot of black regions (threads remaining idle). These inputs from the Paraver visualization helped us to determine some code transformations to make better use of unoccupied computational resources. The key point was to point out the cause of the problem, the improvements were not so difficult to put into place.
CONCLUSION gysela achieves a good parallel computation scalability up to 64k cores combining several levels of parallelism and a hybrid OpenMP/MPI approach. The transpose algorithm helps to reduce the CFL and represents a scalable and robust solution to handle different physics regimes. Global execution times are greatly reduced thanks to this upgrade. Simultaneous multi-threading of Haswell architecture enables a program to better use each core efficiently with multiple threads. Therefore, by using multiple threads on a single core, we are able to increase the cumulated throughput of all available integer units and floating-point units. In gysela code, it leads to a significant reduction of the execution time, larger than 20% in typical cases. Some adaptations of MPI and OpenMP usage to avoid synchronization that conflicts with threads and SMT further improve performance.
This article reports on the development and tests of the adaptive semi-implicit scheme (ASIS) solver for the simulation of atmospheric chemistry. To solve the ordinary differential equation systems associated with the time evolution of the species concentrations, ASIS adopts a one-step linearized implicit scheme with specific treatments of the Jacobian of the chemical fluxes. It conserves mass and has a time-stepping module to control the accuracy of the numerical solution. In idealized box-model simulations, ASIS gives results similar to the higher-order implicit schemes derived from the Rosenbrock's and Gear's methods and requires less computation and run time at the moderate precision required for atmospheric applications. When implemented in the MOCAGE chemical transport model and the Laboratoire de Météorologie Dynamique Mars general circulation model, the ASIS solver performs well and reveals weaknesses and limitations of the original semi-implicit solvers used by these two models.
This solver can be very efficient but it is not naturally mass conserving. It is, for example, implemented in the CHIMERE model  (Menut et al. , 2013). Mass-conserving, multistep or multistage, and high-order accurate implicit methods exist to solve the ODE stiff system. Among the methods based on BDF, Gear's predictorcorrector method has been adapted to atmospheric chemical systems, for example, the SMVGEAR code  (Jacobson and Turco, 1994)  implemented in the GEOSCHEM CTM  (Bey et al. , 2001). More recently, the Rosenbrock's method  (Rosenbrock, 1963)  is becoming widely used in atmospheric chemistry modeling  (Sandu et al. , 1997b)  despite the fact that its computational cost is still rather high compared to approaches based on low-order BDF methods. The implementation in chemical models of Rosenbrock's and other highorder methods has been eased by the development of the Kinetic PreProcessor (KPP) by  Sandu and Sander (2006) , which allows the choice of an integration method and generates the adequate codes accordingly.
An example is given in Sect. 5 with the simulation of the atmospheric composition of the planet Mars. Once the matrix M is evaluated and the time step is determined (see the next section), the solver computes the solution to the system of linearized Eq. (  7 ), which becomes a possible computational bottleneck. Our approach is to use standard methods and well-optimized software libraries. Our baseline option is to use the direct solver DGESV of the Lapack library that solves system (6) by lower-upper (LU) decomposition. Therefore no extra specific routine associated with the chemical mechanism is needed and the optimization on the computer used is left to the implementation of the Lapack library. As reported below, this option works well and gives accurate results even for comprehensive mechanisms involving hundreds of species or more.
The FLUX case To assess the performances of ASIS, two reference simulations have been obtained for the FLUX case using Rosenbrock's and Gear's BDF solvers (referred to hereafter as R1 and G1). Those solvers use the "ode23s" and "ode15s" codes, respectively, from the Matlab ODE suite  (Shampine and Reichelt, 1997; Ashino et al. , 2000). For the Rosenbrock's scheme a three-stage algorithm is used and the simulations are third-order accurate. For the Gear's scheme the third-order accurate option was also chosen. In the most accurate simulations, R1 and G1 (see Table  4 ), the relative tolerance RTOL is set to 0. 001 and the absolute tolerance ATOL equals 10 4 molecules cm -3 for all species. With the Rosenbock's solver additional simulations with higher tolerance values have been performed, R2 with RTOL = 0.
01 and R3 with RTOL = 0. 025. The same FLUX case is integrated using the ASIS solver. In a first simulation noted A1, ASIS uses a RTOL value of 0. 001 and a minimum time step of 1 s. For the solution of the linear system associated with ASIS, the DGESV code of the Lapack library is used. To compare with the Rosenbrock's solver a second simulation A2 has been obtained with the same settings as A1 but with a higher relative tolerance value of 0. 01, and a third one A3 with a tolerance value of 0. 025. For all experiments the FLUX case is integrated over 24 h. The settings used in the overall simulations are given in Table  4. Figures  1  and 2  show the evolution of some key species for each experiment and the relative differences from the R1 experiment. Those results are representative of all of the species.
The efficiency of ASIS in terms of CPU time has been evaluated within the Matlab environment. Table  4  gives the ratio of CPU time used for each simulation relative to the R1 case. The ode23s code used to run the R cases needs the implementation of two subroutines, one that computes the species tendencies and one that gives the Jacobian of the system. If the latter is not provided, the ode23s code computes an approximation of the Jacobian by differentiation and the CPU cost increases by a factor of 2 to 10. As can be seen from Table  4 , the CPU cost of ASIS is comparable to or lower than the ode23s cost for relative tolerance values larger than 0. 01. An important point to mention is that within the Matlab environment the CPU cost does not come from the linear algebra parts of the algorithms but from the evaluation of tendencies and Jacobian matrices.
Therefore it is very dependent upon the chemical system and the details of the programming of the associated subroutines. The situation is quite different within the Fortran environment. With the Fortran version of ASIS the CPU cost for the calculation of the approximated Jacobian (the matrix M of Eq. 7) is negligible compared to the linear algebra computations. This is because the compiler efficiently handles the associated subroutine (fill_matrix, see Sect. 7) that contains frequent indirect addressing. It is not possible to evaluate if this is also the case with all the codes based on Rosenbrock's algorithm, but if it is so, ASIS should perform well when the mean time steps are comparable since it needs fewer linear algebra computations. For the A experiments, ASIS uses the DGESV code for the solution of the linear systems.
Conclusions The ASIS solver has been designed to cope with the various situations encountered within the numerical simulation of the atmospheric chemistry. The main properties of the solver are mass conservation, an approximation of the Jacobian matrix of the chemical fluxes that stabilizes the associated system of differential equations, a time-stepping varying module to control accuracy, and a code implementation that allows an easy adaptation to various chemical schemes. In box-model test cases, the numerical solutions obtained with the ASIS solver were found to be in good agreement with those of multistep algorithms like Rosenbrock's and Gear's methods. The ASIS solver has been implemented in 3-D models of the Earth (MOCAGE) and Mars (LMD model) planets. The results with MOCAGE using ASIS reveals two weaknesses of the original semi-implicit solver.
For some applications the aqueous phase associated with the presence of clouds must be also considered (e. g. ,  Leriche et al. , 2013). The chemistry module has to solve both gaseous-and aqueous-phase chemistry as well as mass transfer reactions between gas and liquid phases. There is a priori no difficulty in adding the prognostic concentrations in the water phase to the system of equations and making a linearization similar to what is done in Eq. (  6 ). However, the addition of aqueous reactions tends to increase the stiffness of the numerical ODE  (Audiffren et al. , 1998) , so the performances of ASIS could diminish and may result in reduced time steps and increased computer time. In conclusion, the ASIS solver can deal with many situations encountered in modeling atmospheric chemistry for a computational cost affordable by CTMs and GCMs that include comprehensive chemical schemes.
Acknowledgements. This work was supported by the  Monitoring Atmospheric Composition and Climate  (http://www. gmes-atmosphere. eu/)and  COPERNICUS Atmospheric Monitoring Service  (https://atmosphere. copernicus. eu/)EU projects. The authors thank  Emanuele Emili ,  Andrea Piacentini , and  Michael Prather  for their very helpful comments on the manuscript. Edited by: G. A. Folberth Reviewed by: three anonymous referees.
In this paper we describe a new rotating module for the PE cell called rotational tomographic Paris Edinburgh cell (RoToPEc), inspired from the stress cell developed in  [7]  and the rotator for single-crystal neutron diffraction at HP detailed in. [8]  The RoToPEc permits controlled torsional shear stresses between the opposed anvils, allowing simultaneous HP and HT measurements of samples subjected to a wide range of strain regimes and provides 360°angular access in the equatorial plane (with a very high accuracy on the angle of anvils rotation) and hence full tomography data sets for in situ HP, HT and/or high stress analyses. Before the RoToPEc, only GSECARS at the Advanced Photon Source at Argonne National Laboratory has demonstrated the capability of performing absorption-based tomography in largevolume samples under HP, HT and high stress using a modified Drickamer cell allowing the full rotation of the sample chamber under load.
[9]  This device has been successfully used to study the volumetric properties of non-crystalline materials,  [10]  and the deformation behaviour of multi-phase composites under shearing. [11]  Nevertheless, this tomography module is installed within a hydraulic press frame on a dedicated beamlime which limits the application to other 3D imaging techniques in various beamlines/synchrotrons, and especially the development of 3D imaging of samples held under extreme conditions on dedicated, high resolution imagining beamlines. In this paper, this new device is fully described and we report the various modifications we have made to adapt the V7 PE press in order to construct our RoToPEc. A detailed description of proof of concept first experiments performed on a high resolution imagery beamline (I12 at Diamond) is then given, which illustrate the efficiency of the set-up and the data quality that can be obtained.
Similar elongation of the PE press columns have been performed previously on the V7 (50%)  [7]  and V4 (88%) models. [13]  Secondly, the piston has been modified to prevent its possible coupled rotation during anvil rotation under load (Figure  1(b) ). In our new system, both anvils can rotate independently under load, with no limitation in angle, through two sets of gear reducers and thrust bearings. A schematic 3D view with cross section of the module is shown in Figure  2. The press including the module has a total weight of 197 kg but the press can be easily disassembled in a few minutes to facilitate transport and reassembly at beamlines (Cf. Figure  3 ), where its integration only requires x-y-z translations and a rotation stage for micro-positioning of the sample along the beam. To illustrate the working principle of the RoToPEc module, some magnified schematic cross-sections of the apparatus are presented in Figures  4  and 5.
Electronics and software control An electronic box has been built following the automation structure detailed in Figure  8. It allows monitoring and control of two motors (FAULHABER 4490 series) rotating both the anvils. It also ensures the reading of the two incremental encoders (Numerik Jena) coupled to the anvil, the encoders motors and port enabling communication between the control software and the electronic box. The PC-interface communication is provided by an RS232 serial link or a Transmission Control Protocol/Internet Protocol Ethernet link and can be installed easily on any computer or laptop. Our home-made user-friendly application on Labview software can operate individually or simultaneously on both motors and send instructions for rotation to specified angles or speed. It also displays continuously the real position of each anvil offering off-line control for the RoToPEc module.