The large processing requirements of seismic wave propagation simulations make High Performance Computing (HPC) architectures a natural choice for their execution. However, to keep both the current pace of performance improvements and the power consumption under a strict power budget, HPC systems must be more energy e cient than ever. As a response to this need, energye cient and low-power processors began to make their way into the market. In this paper we employ a novel low-power processor, the MPPA-256 manycore, to perform seismic wave propagation simulations. It has 256 cores connected by a NoC, no cache-coherence and only a limited amount of on-chip memory. We describe how its particular architectural characteristics influenced our solution for an energy-e cient implementation. As a counterpoint to the low-power MPPA-256 architecture, we employ Xeon Phi, a performance-centric manycore. Although both processors share some architectural similarities, the challenges to implement an e cient seismic wave propagation kernel on these platforms are very di↵erent. In this work we compare the performance and energy e ciency of our implementations for these processors to proven and optimized solutions for other hardware platforms such as general-purpose processors and a GPU. Our experimental results show that MPPA-256 has the best energy e ciency, consuming at least 77 % less energy than the other evaluated platforms, whereas the performance of our solution for the Xeon Phi is on par with a state-of-the-art solution for GPUs. Public policies for risk mitigation and damage assessment in hypothetical 2 earthquakes scenarios as well as oil and gas exploration make an intensive use 3 of simulations of large-scale seismic wave propagation. In order to provide re-4 alistic results, these simulations rely on complex models, which often demand 5 considerable computational resources to reliably process vast amounts of data 6 in a reasonable amount of time. Such processing capabilities are often provided by High Performance Com-8 puting (HPC) platforms. Up until recently, HPC platform capabilities were 9 evaluated almost exclusively based on their raw processing speed. However, one 10 of the aspects that hinder the quest for ever growing performance is the exces-11 sive use of energy. For that reason, the energy e ciency of HPC platforms has 12 become, in some contexts, as important as their raw performance. Indeed, the seek for alternatives to lower current energy consumption arose 14 first in the embedded systems community but lately became one of the ma-15 jor concerns of the HPC scientific community [1, 2]. Recently, manycore pro-16 cessors, a new class of highly-parallel chips, was unveiled. Tilera Tile-Gx [3], 17 Kalray MPPA-256 [4], Adapteva Epiphany-IV [5], Intel Single-Chip Cloud Com-18 puter (SCC) [6] and Xeon Phi are examples of such processors, providing up to 19 hundreds of autonomous cores that can be used to accomplish both data and 20 task parallelism. This distinctive characteristic sets them apart from SIMD-like 21 highly-parallel architectures such as Graphics Processing Units (GPUs). For 22 this reason, in this paper, we classify GPUs separately, in their own category. While some manycore processors may present better energy e ciency than 24 state-of-the-art general-purpose multicore processors [7], their particular archi-25 tectural characteristics make the development of e cient scientific parallel ap-26 plications a challenging task [5, 8]. Some of these processors are built and opti-27 mized for certain classes of embedded applications like signal processing, video 28 decoding and routing. Additionally, processors such as MPPA-256 have impor-29 tant memory constraints, e.g., limited amount of directly addressable memory 30 and absence of cache coherence protocols. Furthermore, communication costs 31 between cores are not uniform, they depend on the location of the commu-32 nicating cores. One of the main reasons for this di↵erence in communication 33 cost is the Network-on-Chip (NoC). When applications distribute work among 34 the cores of the processor, they should take into consideration the placement 35 of these tasks and their communication patterns to minimize communication 36 costs. This problem is similar to that faced by applications running on NUMA 37 platforms [9, 10], or MPI applications that try to take into account both the 38 network topology and the memory hierarchy to improve performance [11, 12]. In this paper we outline the architectural distinctiveness of the low-power 40 MPPA-256 manycore processor. Considering these characteristics, we describe 41 how the main kernel of a seismic wave propagation simulator was adapted to 42 h.aochi@brgm.fr (Hideo Aochi), navaux@inf.ufrgs.br (Philippe O. A. Navaux), jean-francois.mehaut@imag.fr (Jean-François Méhaut) 2 this platform. Due to the limited size of the local memories in MPPA-256, we 43 developed a new multi-level tiling strategy and a prefetching mechanism to 44 allow us to deal with real simulation scenarios and to alleviate communication 45 overheads. We also describe the di culties and solutions (some of them generic 46 enough to be used in di↵erent contexts) we employed during this adaptation. As a counterpoint to the low-power MPPA-256, we also adapted the same 48 kernel to a performance-centric manycore, the Xeon Phi. Although both pro-49 cessors share some architectural characteristics such as the presence of a NoC, 50 non-uniform costs of communication between cores and limited fast local mem-51 ories, they have their own peculiarities which directed us to employ di↵erent 52 approaches in each case. For instance, in both architectures, each thread has 53 access to approximately 128 kB of fast local memory. However, their memory 54 organization is very di↵erent. On Xeon Phi each core has 512 kB of L2 shared 55 by up to 4 threads. On the other hand, on MPPA-256 the processing cores have 56 8 kB of L1 cache each and cores are grouped into clusters of 16. Each one of 57 these clusters has 2 MB, of work memory shared among the cores. Contrary to 58 Xeon Phi, MPPA-256 does not have cache-coherence or automatic load/store of 59 values from/to main memory, and all memory transferences must be handled 60 by the programmer explicitly. Ultimately, this means that in order to be able 61 to achieve good performance on MPPA-256, the application is obliged to tightly 62 control communications between the cores in the same cluster, between clusters, 63 and between each cluster and the main memory. This paper explains the de-64 tails that made us decide for a multi-level tiling strategy proposed for MPPA-256 65 instead of a more traditional approach of cache blocking that was used for Xeon 66 Phi. 67 Taking as a basis optimized GPU and general-purpose processor implemen-68 tations for this kernel, we show that even if MPPA-256 presents an increased 69 software development complexity, it can indeed be used as an energy-e cient 70 alternative to perform seismic wave propagation simulations. Our results show 71 that the solution we propose for the MPPA-256 processor achieves the best en-72 ergy e ciency, consuming 77 %, 86 % and 88 % less energy than optimized so-73 lutions for GPU, general-purpose and Xeon Phi architectures, respectively. The 74 performance achieved on Xeon Phi, on the other hand, is comparable to a state-75 of-the-art GPU. Moreover, the execution was 58 % and 73 % faster than that of 76 multicores and MPPA, respectively. 77 The remainder of this paper is organized as follows. Section 2 discusses 78 the fundamentals of seismic wave propagation. Then, Section 3 discusses the 79 main challenges we must overcome when dealing with parallel seismic wave 80 propagations on MPPA-256 and Xeon Phi. Moreover, it presents our approaches 81 to perform e cient seismic wave propagation simulations on both processors. 82 Section 4 discusses performance and energy e ciency experimental results. Sec-83 tion 5 describes related work and Section 6 concludes this paper. 84 2. Seismic Wave Propagation 85 Quantitative earthquake hazard assessment is crucial for public policy, re-86 duction of future damages and urban planning. Recent important geological 87 3 events, such as the M W 7.9 (moment magnitude) earthquake in China (2008) or 88 the M W 8.8 earthquake in Chile (2010), have been studied using numerical tools 89 with good agreement with observation. From a physical point a view, the mech-90 anisms used to describe earthquake events are complex and highly nonlinear. 91 For the sake of simplicity, wave propagation simulations usually consider elastic 92 domains and include some anelasticity to account for attenuations. However, 93 in the case of strong motion, these simplifying assumptions are not capable of 94 correctly describing the phenomenon. A fundamental challenge for earthquake 95 engineering is predicting the level and variability of strong ground motion from 96 future earthquakes. Enhancing this predictive ability requires understanding of 97 the earthquake source, the e↵ects of the propagation path on the seismic waves, 98 and basin and near-surface site e↵ects. In this paper, we restrict our discussion to the elastodynamics equations 100 corresponding to the elastic situation using a linear model. However, our ap-101 proach could also be used in nonlinear models that rely on a linear solution 102 inside the iterative procedure such as Newton-like methods. Among numeri-103 cal methods, variational methods such as finite element method [14] or spectral 104 element method [15] can be employed for both regional and global scale seismol-105 ogy whereas the Discontinous Galerkin method is very e cient for earthquake 106 source modeling [16]. In terms of popularity, the Finite Di↵erence Method 107 (FDM) remains one of the top contenders due to its implementation simplicity 108 and e ciency for a wide range of applications. Although this method was ini-109 tially proposed almost thirty years ago [17], this classic discretization technique 110 has been the focus of continuous refinements since its proposal [18, 19]. 111 Simulations of seismic wave propagations are often constrained by the com-112 putational and storage capacity of the hardware platform. Thus seismologists 113 often reduce the scope of the simulation to a small volume instead of considering 114 the whole Earth. The memory footprint of each simulation depends on both 115 the domain of interest and the number of grid points per wavelength. This last 116 condition guarantees numerical stability and therefore the accuracy of the sim-117 ulation. In our case, for practical reasons, we limit our simulations to problem 118 domains that fit into the 2 GB of memory available on the MPPA-256 platform. 119 14 e cient utilization of the caches through cache blocking alone increased the 479 performance by 79 % on average. If used in isolation, cache blocking does not 480 allow a full use of Xeon Phi's vectorization capabilities and creates additional 481 cache faults due to unaligned memory accesses. By using these two techniques 482 at once we can, at the same time, make a good utilization of Xeon Phi's vectorial 483 capabilities and its fast local caches. Indeed, when these two techniques are 484 used together we see an average improvement of 92 % in performance, and an 485 improvement of 91 % for 224 threads. 486 Dynamic and guided scheduling policies had a better scaling behavior as well 487 as better performance than static scheduling. On average, dynamic scheduling 488 provided an increase of 14 % over static scheduling. Static scheduling has the 489 downside of causing imbalances during the execution of the simulation with an 490 intermediate number of threads. For instance, for 58 threads, iterations would 491 be divided equally between threads, but since a single core would execute two 492 threads, it would actually end up executing twice as many iterations as the other 493 cores. This behavior can clearly be seen as steps in the speedup graph. However, 494 this assumption is not enough to e↵ectively explain why the performance using 495 static scheduling is not on par with (or better than) the performance of dynamic 496 scheduling when the number of threads is a multiple of the number of cores. The reason for that lies on the high usage levels of the NoC during the 498 implicit synchronization between the calculation of the velocity and stress com-499 ponents of the stencil. For each timestep, the stress calculation loop only starts 500 its execution after all the velocity components of the previous loop have finished 501 (cf. Algorithm 2.1). This synchronization makes every thread to try to commu-502 nicate at the same time at the start of each parallel loop. On the static scheduler, 503 the concurrency for the NoC might introduce small delays for the beginning of 504 the execution of each thread, essentially making the loop's critical path slightly 505 longer. While by itself each one of these delays is not very significant, they 506 accumulate at each timestep both at the stress and velocity calculation loops. Since the concurrency for the NoC also increases with the number of threads, 508 the relative execution performance di↵erence between the static and dynamic 509 approaches also increases. This imbalance is compensated on the dynamic case, 510 even if the dynamic scheduler has a higher overhead, because it is able to keep 511 more threads active up to the end of each loop. Guided scheduling also allows for a better work distribution and smaller 513 imbalances than the static scheduler between threads at the end of each loop. However, the synchronization between velocity and stress loops and the small 515 delay introduced during the initial phases of each loop still create slight imbal-516 ances at the end of both the stress and velocity loops of each timestep. Using 517 the guided scheduler we were able to achieve near linear speedups up to the 518 number of cores. After this point, the dynamic scheduler gradually approaches 519 the guided performance to eventually beat it with the best execution time with 520 224 threads. Finally, for static scheduling, thread a nity was able to improve the perfor-522 mance of MA+CB in 11 %, on average. On the other hand, by its nature, dy-523 namic and guided scheduling reduce the importance of thread a nity. Indeed, 524 interesting than The Altix UV 2000 platform has the best trade-o↵ between the time-to-596 solution and energy-to-solution, for both maximum performance and optimal 597 energy consumption settings. If we consider ED 2 P, this di↵erence becomes 598 even more noticeable. The significant di↵erences highlighted in the combined 599 time and energy evaluation can be explained by the poor arithmetic intensity 600 Using the low-power Adapteva Epiphany-IV manycore, Varghese et al. [5] 633 described how a stencil-based solution to the anisotropic heat equation using a 634 two-dimensional grid was developed. This manycore has a low power budged 635 (2 W) and has 64 processing cores. Each core has a local scratchpad of 32 kB and 636 no cache-memory. Cores are connected by a NoC which allows direct memory 637 access to 1 GB of RAM and also to other cores' local memories. The challenges 638 we faced for the development of our seismic wave propagation simulation on the 639 MPPA-256 architecture are similar to the ones faced by the authors of this work. 640 Similar to MPPA-256, Epiphany-IV has a very limited amount of local memory 641 available to each core and no automatic prefetching mechanism exists; every 642 data movement has to be explicitly controlled by the application. The authors 643 demonstrated that, even if the implementation of an e cient application can 644 be a challenging task, their solution was able to achieve a FLOPS/Watt ratio 3 645 times better than an Intel 80-core Terascale Processor. 646 Adapting seismic wave propagation numerical kernels to emerging parallel 647 architectures is an active research topic. Due the versatility of the FDM used 648 both for oil and gas applications and standard earthquakes modeling, several 649 strategies have been proposed. For instance, Intel and NVIDIA have proposed 650 optimized implementations of such stencils [25, 32]. 651 Dupros et al. [33] presented a review of the scalability issues for distributed 652 and shared-memory platforms with a focus on mapping processes/threads on 653 hierarchical clusters. Dursun et al. [34] introduced several additional strategies, 654 including inter-node and intra-node optimizations. Recently, Christen et al. [35] 655 described the use of the Patus framework to optimize the AWP-ODC finite 656 di↵erence code. In particular, they detailed the impact of vectorization and 657 cache prefetching and reported a performance improvement up to 15 % when 658 running the complete application with 512 MPI processes. 659 Several research e↵orts have been done on the adaptation of seismic wave ker-660 nels on GPUs to overcome the poor arithmetic intensity (FLOPS/transferred 661 byte ratio) o↵ered by elastodynamics stencils [27, 36]. The results obtained 662 demonstrate that GPUs could be a relevant alternative to standard general-663 purpose processors. Krueger et al. [37] compared the performance and energy 664 e ciency of an Intel Nehalem platform, an NVIDIA Tesla GPU and a simulated 665 general-purpose manycore chip design optimized for high-order wave equations 666 called "Green Wave". They showed that Green Wave can be up to 8x and 3.5x 667 more energy e cient per node when compared with the Nehalem and GPU 668 platforms, respectively. Di↵erently from this work, our experiments and mea-669 surements were carried out on real manycore processors. Even if the low-power 670 MPPA-256 is not optimized for high-order wave equations as Green Wave, we 671 were still able to achieve energy e ciency improvements of at least 86 % and 672 77 % when compared to the Xeon E5 and Tesla K20, respectively. 673 Conclusion and Perspectives The use of highly-parallel manycore processors to meet the high demand for 675 data processing capabilities needed by parallel scientific simulations has recently 676 19 become one of central points of interest in the HPC community. In this paper, 677 we presented our approach to the seismic wave propagation simulation using the 678 MPPA-256 and Xeon Phi manycore processors. Although both processors share 679 some architectural characteristics, they present di↵erent challenges that must 680 be overcome in order to obtain good performance results.
Research efforts such as the presented in Castro et al. REF improved and evaluated the performance of the acoustic wave propagation equation on Intel Xeon Phi and compared it with MPPA-256, general-purpose processors and a GPU.