A strategy is proposed for the development of algorithms and software for Exaflops supercomputers. The strategy contains three stages. The first stage is the co-design, which is defined as considering the architecture of the supercomputer at all stages of the development of the code. The second is the forward-looking development of algorithms and software for the most promising Exaflops supercomputers. The forward-looking development is based on the simulation of the algorithm behavior within the given supercomputer architecture. The third stage is the estimation of the energy efficiency of the algorithm with different implementation for one single architecture or with different supercomputer architectures. The strategy is shown by the example of two problems from astrophysics and plasma physics.
The bottleneck in the use of the Exaflops supercomputers will be the incredible degree of parallelism, hundreds of millions or even billions of simultaneously working computation cores. Software development for supercomputers is closely connected with the evolution of their architectures. Thus the study of scalability of the parallel algorithms with the future Exaflops supercomputers is an urgent question. One must note that the computational algorithms evolve slowly compared to the evolution of hardware. So, even for Exaflops supercomputers the present-day computational algorithms will be applied since there is no special mathematics that takes the peculiarities of Exaflops machines into account. Therefore it is important to study right now the perspective of the implementation of the present-day algorithms on future supercomputers and to investigate their parallelism, scalability and energy consumption.
In [
The problem of simulation of the algorithms scalability is not something new; there is a lot of
work throughout the world [
In the Institute for System Programming RAS a model of a parallel program was developed under
the supervision of Academician V.P.Ivannikov. The model can be efficiently interpreted by an
instrumental computer thus providing the precise prediction of the real execution of a parallel program
within the given computer. The model is developed for the parallel programs with explicit message
interchange using Java together with MPI. Therefore the model is included into ParJava [
In [
Since the architecture of the exaflops supercomputer is an open question, the energy efficiency of
the algorithms becomes a question of great importance. Energy efficiency may be dividing into 3 main
parts [
Co-design of the parallel methods for the solution of large-scale problems is a process complex enough and hard to formalize. This is because every problem has its own features. Even if two problems are from the same area of science, they will still have some differences. The careful consideration of these features and differences may alter not only the parallel computational method, but the mathematical model itself and of course, the program implementation and its efficiency.
It is impossible, of course, to make a «collection of recipes» for efficient solution of all the largescale problems. However, some general approach is possible to be proposed. In any case, the concept of co-design consists of the following: 1. The formulation of the physical statement of the problem. 2. The mathematical formulation of the physical problem.
3. The development of the numerical method for the solution of the mathematical formulation of the physical problem.
4. The selection of data structures and of the parallel algorithm.
5. The consideration of supercomputer architecture.
6. The use of the development tools.
At the stage of the formulation of the physical statement of the problem it is possible to define the
main physical process. This will be the basis on the following stages. For example, the simulation of
gas flows may be reduced to the solution of hyperbolic equations (gas dynamics equations) in the case
of jet pipe simulation [
Even at the stage of the formulation of the physical statement of the problem it is possible to
begin using some parallel formalism and to define the type of parallelism that will use in the solution
of the problem. To our opinion the most natural division of the types of parallelism is the following:
distributed problems and the problems with dependence parallelism. With no loss of generality all the
distributed problems are integration problems. It means that the problem is decomposed into the set of
independent subproblems, the results of which are accumulated in some shared object. In the case of
numerical integration or summation in quantum chemistry problems [
The problems with dependence parallelism may be divided into three types: explicit operator
inversion problem; implicit operator inversion problems solved by grid method; implicit operator
inversion problems solved by grid-less methods. Operator inversion problems are primarily the solution of
either elliptic equations [
Implicit operator inversion by the grid methods is the solution of hyperbolic equations within a
compact stencil for most of the continuous medium models: gas dynamics equations [
The present authors develop the Lagrangian-Eulerian approach for solving astrophysical problems on the basis of Fluids-In-Cells (FlIC) method and Godunov method for several years. Using of this combination enables to eliminate the Galilean non-invariance and the use of Godunov method at the Eulerian stage enables to simulate the discontinuities correctly.
GPU architecture is a set of concurrently executing threads, combined in a multi-level 2d topology. The maximum performance on this architecture can be obtained in the absence of synchronization between threads. In this case the strategy to co-design a numerical method is to use the computational mesh topology, which directly projected on the topology of GPU architecture (e.g. use of regular mesh), and independent computation of the values at the next step in each cell of the computational domain. The Riemann problem for Euler stage original method satisfies this property, but for the Lagrangian stage needed a modification. This modification allows independently computing advective transport for each cell of the computational domain.
Are traditionally, collisionless component is described using N-body models. Nevertheless, this model has disadvantages as a spurious generation of entropy, increased communication overhead and poor load balancing. Therefore, in co-design algorithms to describe collisionless components in the context of modeling the galaxies collision was selected «Collisionless Stellar Hydrodynamics» approach. This approach is based on the equations for the first moments of the collisionless Boltzmann equation. This approach is adequate for importance of a movement particles cluster and high kinetic energy. In addition, for describing «Collisionless Stellar Hydrodynamics» models can be used an unified numerical methods and parallel algorithms, that are used to solve the hydrodynamic equations. Consequently, they can be efficiently implemented on GPU-based supercomputers.
In the present case co-design begins at the stage of the physical consideration of the problem. It is known from experiment that the plasma density modulation cannot exceed 300 \%. It means that the number of model particles cannot be increased without a limit. Thus there is no need for dynamic load balancing, and the absence of dynamic balancing improves the reliability and efficiency of the parallel implementation.
At the stage of the numerical method design the field evaluation method was chosen that is built on the basis of Faraday and Ampere laws. In such a way there is no need to solve Poisson equation. Instead, the equations that represent the Faraday and Ampere laws in the numerical form are solved by the Langdon-Lasinski scheme. This results in the field solver with virtually unlimited scalability.
At the stage of supercomputer architecture selection the PIC method details are taken into account. In order to evaluate the new values of position and impulse of a particle it is necessary to know the values of electric and magnetic fields at the present position of the particle. Each of the three components of both electric and magnetic field is stored in a separate 3D array. In such a way six 3D arrays are accessed at each time step for each particle. Since the particles are situated randomly in the domain, the access to the arrays is also unordered. It means that the use of the cache memory can not reduce the computation time. If a part of the field array was fetched to the cache in the process of computation of the particle movement, it would be impossible to use this part of the field array for the computation with the next particle, because it is (most likely) situated in a completely different part of the subdomain. Since the cache memory cannot store all the six arrays for fields, one has to access the RAM (Random Access Memory) for computation of the particle movement. And since the performance of the processor is usually limited by the memory bandwidth, it is the memory bandwidth that determines the speed of the computation with particles and the performance of the program as a whole (particles take from 60% to 90% of the total time).
This fact determines the transition to the supercomputers equipped with GPUs because CUDA has a lot of tools to accelerate memory use for the PIC method.
At the stage selection of software design tools is the following. For the PIC method with a very big number of independently processed elements (the model particles) the use of CUDA technology is very efficient. Other parallel technologies for hybrid supercomputers such as OpenCL, OpenMP, OpenACC could be also used but it is CUDA that gives the possibility to use the highest number of parallel processes and to get the highest performance.
The last stage of co-design is the adaptation of the algorithm to the GPU architecture. The traditional PIC method implementation — all the particles stored in one very big array is unacceptable for GPUs because it prevents the use of the advantages of GPU memory. Thus the particles are distributed between the cells and the computation is conducted as follows: one block of threads for one cell.
The multi-agent simulation system AGNES [
AstroGrid Agent sending messages to two neighbors and expects to receive the messages from both neighbors as implemented in software. A calculation continues only after this cycle.
PlasmaGrid agents are summarizing data for all agents before the next iteration of the loop.
All data sending delays are collected and analyzed according to real run tests of astrophysics and plasma physics codes.
The model of interacting processes is defined here as the set of threads (not Operating System threads), being concurrently executed on a single computational node. These threads interact with each other by means of information interchange, for example, with MPI. The main parameter of the thread is the execution time and also the time spent on passing the data to another thread. Considered this, let us give the scheme of interactions between threads for the problems of astrophysics and plasma physics.
The peculiarity of parallel grid methods is the possibility for geometrical decomposition of computational domain and further exchange by the boundary values by the adjacent computational nodes only. These methods are the main tool for the solution of hyperbolic equations, not only in astrophysics. Here, if the computational domain contains N 3 cells, then the number of element transmitted to the adjacent computational node is N 2 elements. The interaction scheme may be shown by the following figure:
stage for N3 cells
N2 cells
stage for N3 cells
N2 cells
stage for N3 cells
N2 cells
stage for N3 cells
N2 cells
stage for N3 cells
N2 cells
stage for N3 cells
N2 cells
In order to investigate the scalability of grid methods with CPUs, GPUs and Intel Xeon Phi accelerators let us consider the following:
1. Execution time for Eulerian and Lagrangian stage of the Fluids-In-Cells method are given by the formulas:
T Euler N 3
E
T Lagrange N 3
L
It means that for each type of the computational device the optimal computation time is known for each stage. The total computation time is proportional to the number of cells.
2. The scalability and speedup within a single CPU or accelerator will not be simulated. The reason is that the number of cores of CPUs and accelerators is fixed and also, this number is much less than the prospective number of computational nodes in peta- or exaflops supercomputers. Therefore the execution time for the stages E,L is optimal for: • 12 cores of Intel Xeon E5-2697 processor • 1024 cores of the Nvidia Kepler GPU • 60 cores of the Intel Xeon Phi 7110 accelerator in the offload mode • 240 cores of the Intel Xeon Phi 7120 in the native mode
In the case of the CPU the speedup value was 12 within one Intel Xeon (12 cores compared to 1 core), for the GPU the speedup was 55 (1024 cores of Kepler compared to 1 core) and for the Intel Xeon Phi accelerator the speedup was 27 in offload mode (60 cores compared to 1 core) and 54 in native mode (240 core compared to 1 core). This speedup is because of the architecture of the task manager of the computational device. The speedup data was obtained by the computational experiments with GPUPEGAS and AstroPhi codes.
3. The communication time is considered to be a linear function of the number of the transmitted element with the latency taken into account:
T Comm C N 2 where is the latency and C is the average time for transmitting of one data element. 4. The feature of the numerical method is that the number of the elements being transmitted is the same for both Eulerian and Lagrangian stages. Consequently, the communication time is also equal.
5. The evaluation and comparison of the communication times is impossible since the computational experiments were conducted with different hardware and different architectures, and in some cases it was just one node of or even one single GPU. Thus it will be presumed that the network of the Siberian Supercomputer Center (SSCC) is being used.
6. In the general case for the solution of the 3D problems with the large number of computational nodes the most efficient geometrical decomposition is the multi-dimensional (among others 3D) decomposition. Most likely, the multidimensional Cartesian topology will be used in Exaflops supercomputers. In figure 1 the linear topology is shown that corresponds the real SSCC architecture. The use of such topology for the interaction between the nodes simplifies the simulation of the execution of threads At the same time it does not limit the use of multidimensional topologies. This is because the data exchange between threads will grow proportionally to the number of dimensions. The communication time for the K-dimensional decomposition will be connected to the communication time along one dimension according to the formula TCDomm T1CDomm T Comm which does not break the generality of the model presented by figure 1.
7. It must be noticed that in the present model only the extensivity of the computer system is being considered. It is quite natural that both computation performance and data exchange speed will grow with Petaflops computing becoming more common and while moving towards Exaflops. But along with this the computational complexity will also be growing, since the physical model of the astrophysical processes will become more complex, resulting in the increasing number of operations per one cell. Due to this reason we investigate the scalability of concurrently executed processes considering the current hardware parameters and consequently, current level of complexity of the computational models. At present, the characteristic number of cell per one accelerator is N 2563 16106 .
In order to describe the model of interacting processes let us write the table of the computation and communication times per one element (one cell).
The goal of this section is to show the parallelization strategy that must be supported by the PIC implementation template since it is quite clear that the plasma physics problems of interest will never fit into a single CPU of GPU memory.
The program was parallelized by the domain decomposition method. The computational domain is divided into parts along the direction orthogonal to the direction of the beam (along the Y axis, the beam moving along the X axis). The computational grid in the whole domain is divided into equal parts (sub-domains) along the Y axis. Each sub-domain is assigned to a group of processors (in the case of a multi-core system a single core would be called a processor, since no hybrid parallelization like MPI+OpenMP is employed, just mere MPI).
Furthermore, the particles of each sub-domain are distributed uniformly between processors of the group with no regard to their position, as it is shown in Figure 2.
The computational domain is divided into 4 subdomains. The particles of each subdomain are distributed between four processors uniformly with no regard to their position. Different symbols (circle, square, diamond, star) denote particles belonging to different processors in the same sub-domain.
Every processor in the group solves the Maxwell equations in the whole sub-domain, and exchanges boundary values of the fields with processors assigned to adjacent sub-domains. Then the movement equations for the particles are solved, and the 3D matrix of the current density and the charge density are evaluated by each processor. But since the processor has only a part of the particles located inside the sub-domain, it is necessary to sum the matrices through all the processors of the group to obtain the whole current density matrix in the sub-domain.
Interprocessor data exchange is performed by the MPI subroutines. The overall communication scheme generally corresponds the scheme given in figure 1 for the astrophysical code. There is only one important difference. The plasma physics code at present has no peer-to-peer communications, just the summation along the whole domain performed by MPI_Allreduce.
The AGNES simulations of astrophysics code was done for computational nodes with Intel Xeon E5-2697 (5632 nodes with 67584 cores), Nvidia Kepler K40 GPUs (1024 nodes with 1048576 cores) and with Intel Xeon Phi 7110 (4096 nodes with 983040 cores). The AGNES simulations of plasma physics code was done for computational nodes with Intel Xeon E5-2697 (3072 nodes with 36864 cores), Nvidia Tesla2090M GPUs (1536 nodes with 786432 cores) and with Nvidia Kepler K40 GPUs (1536 nodes with 4.4 million cores).
The results for astrophysical code simulation are the following: 1) Running time increases insignificantly (about 20%) up to the 5120 computing nodes; 2) The best scalability was demonstrated for nodes with Nvidia Kepler K40 and Intel Xeon Phi (native mode) accelerators.
The results for plasma physics code simulation are the following: 1) The substantial increase of communication messages gives the 58% growth of code execution time for 3072 nodes with Intel Xeon CPUs;
2) The plasma physics code has the best scalability on Nvidia Kepler K40 GPUs nodes;
In order to improve the energy efficiency of the algorithms for the problems of astrophysics and plasma physics each code was being optimized along the following three directions: the efficiency of interaction with the CPU and the RAM, the efficiency of working with network devices and also the efficiency of working with peripherals, for example, the input-output system. The interaction with the CPU and RAM was facilitated with Intel Parallel Studio XE 2016 Beta package, including Intel Vectorization Advisor. The latter analyzes the code and gives offers for code optimization for maximally efficient use of vector instructions. The use of this package gave the 30% performance gain for the astrophysical code for Intel Sandy Bridge CPU. For both astrophysics and plasma physics problems presented here the most efficient methods were selected. Efficiency here means the minimal number of network data exchanges.
The optimization of the numerical method of the astrophysical code and the use of the last version of the Intel Parallel Studio XE 2016 Beta package resulted in the excellent result for the load balancing for the processor cores. The time for the MPI operations does not exceed 2-3% of the total execution time. Figures 5-7 present the profiles of the codes.
The figures 7-9 show the profile for plasma physics code. The co-design of the numerical algorithm and the use of latest CUDA 7.0 optimization techniques together with Intel Parallel Studio resulted in good load balancing with MPI operation time not exceeding 10%.
In this paper, we proposed a methodology for the development of algorithms and software for exaflops supercomputers. Our approach has three related stages: co-design stage, scalability simulation, energy efficiency improvements. The co-design technique is based on the idea that algorithm and software design takes into account supercomputer's architecture. One of the most important stages in our approach is the scalability simulation. We use AGNES software package to predict the scalability of algorithms on possible exaflops supercomputer's architecture. Energy efficiency improvements stage consists of different code optimizations such as amplification of CPU or GPU usage, load balancing and MPI operations reduction. Our approach is illustrated by two examples: astrophysics and plasma physics codes. As a result of our approach for software development we got highly parallelized codes with suitable load balancing and MPI operations time.