Plasma-PEPSC D5.2 Report on Algorithmic Changes and Memory Optimisation for Heterogeneous Systems in Plasma-PEPSC Codes

HORIZON JU Research and Innovation Actions HORIZON-EUROHPC-JU-2021-COE-01-01 European High-Performance Computing Joint Undertaking Plasma Exascale-Performance Simulations CoE 101093261 D5.2 Report on Algorithmic Changes and Memory Optimisation for Heterogeneous Systems in Plasma-PEPSC Codes WP5: Accelerated Plasma Simulations on Heterogeneous Systems Date of preparation (latest version): 31/12/2024 Copyright©2023 – 2027 The Plasma-PEPSC Consortium D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 2 DOCUMENT INFORMATION Deliverable Number D5.2 Deliverable Name Report on Algorithmic Changes and Memory Optimisation for Heterogeneous Systems in Plasma-PEPSC Codes Due Date 31/12/2024 Deliverable lead HZDR Authors Michael Bussmann (HZDR), Jeremy Williams (KTH), Markus Battarbee (UoH), Etienne Renault (SiPearl) Leon Kos (UL), Jordy Trilaksono (MPG), Gabin Schieffer (KTH) and Felix Jung (TUM) Responsible Author Michael Bussmann (HZDR) E-mail: [email protected] Keywords GPU Accelerators Heterogeneous Memory Systems WP/Task WP5/Task D5.2 Nature R Dissemination Level PU Final Version Date 31/12/2024 Reviewed by Celine Scetbun (SiPearl) Aleˇs Podoln´ık (IPP CAS) and Jakub Hromadka (IPP CAS) D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 3 DOCUMENT HISTORY Partner Date Comment Version KTH 03/03/2024 Skeleton version of the deliverable 0.1 HZDR 02/12/2024 First draft 0.2 KTH 03/12/2024 Final draft updated for internal review 0.3 HZDR 13/12/2024 Revised draft after internal review 0.4 KTH 20/12/2024 Final cleanup for submission 1.0 D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 4 Executive Summary This document provides insight into the current state of porting efforts for the PlasmaPEPSC codes toward using accelerator hardware such as graphics processing units (GPUs). It summarizes state-of-the-art techniques for accelerator programming, porting strategies from CPU to accelerators and usage of multiple GPUs. One focus is the use of heterogeneous memory when using accelerators, showcasing some of the libraries provided by Plasma-PEPSC members for memory management on heterogeneous devices. Although this report documents work in progress, it is evident that almost all PlasmaPEPSC codes are ready for use with GPU accelerators in production on EuroHPC systems. We thus are confident to be able to present regular updates on GPU performance within the project. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 5 Contents 1 Introduction 6 2 Background on Accelerators and Heterogeneous Memory 6 3 Porting Status and Optimization Potential of the Codes 7 3.1 BIT ...................................... 7 3.1.1 Hybrid MPI and OpenMP/OpenACC BIT1 . . . . . . . . . . . . 9 3.1.2 Accelerating BIT1 Asynchronously on Vega and MN5 . . . . . . . 10 3.2 GENE/GENE-X ............................... 13 3.2.1 GENE................................. 13 3.2.2 GENE-X................................ 13 3.3 PIConGPU .................................. 18 3.4 Vlasiator.................................... 19 4 GPU Acceleration Strategies 20 4.1 Programming Interfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.1.1 Multi-GPU Communication Interfaces . . . . . . . . . . . . . . . 21 4.2 Performance Optimisation Strategies . . . . . . . . . . . . . . . . . . . . 24 4.2.1 GPU Matrix Processing Units . . . . . . . . . . . . . . . . . . . . 25 4.2.2 Mixed Precision on GPUs . . . . . . . . . . . . . . . . . . . . . . 27 5 Heterogeneous Memory Approaches 28 5.1 Heterogeneous Memory Systems . . . . . . . . . . . . . . . . . . . . . . . 28 5.1.1 Device Memory and Interconnect Study . . . . . . . . . . . . . . 28 5.1.2 Memory exposure and performance of HBM/DDR . . . . . . . . . 31 5.1.3 GPUDirectStorage.......................... 35 5.2 Programming Interfaces for in-memory data Management in HMS . . . . 36 5.2.1 Concurrent memory allocation with mallocMC . . . . . . . . . . . 36 5.2.2 Optimum memory layout and reduced representations with Llama 38 5.2.3 H2M.................................. 39 5.2.4 StarPU ................................ 39 5.2.5 Umpire ................................ 39 6 Conclusion 39 D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 6 1 Introduction This report provides an in-depth examination of the adaptation of Plasma-PEPSC codes to accelerators and heterogeneous memory systems, with a focus on utilizing advanced EuroHPC high-performance computing (HPC) resources. By leveraging cutting-edge GPU technology and complex memory hierarchies, the report addresses the challenges and opportunities associated with achieving scalable and efficient plasma simulations. The report begins by offering a comprehensive background on the evolution and current state of accelerator hardware, programming models, and heterogeneous memory systems. It explores the transformation of GPUs from specialized graphics processors to powerful tools for scientific simulations. The discussion highlights the critical need for adapting and the simulation codes to these cutting-edge technologies to optimize performance and scalability on modern HPC platforms. It then delves into the current state of implementation, examining the progress made in porting and optimizing key Plasma-PEPSC codes. Detailed analyses of individual codes provide insights into their specific advancements, challenges, and areas for future improvement, showcasing the collaborative effort to overcome technical bottlenecks. An exploration of GPU acceleration strategies follows, presenting the methods and tools employed to harness the full potential of GPUs. This includes a discussion of programming interfaces, multi-GPU communication techniques, and lessons learned in optimizing data movement and computation for heterogeneous systems. The section emphasizes the importance of efficient resource utilization to achieve high-performance outcomes. The report also investigates performance optimization strategies, focusing on critical path analysis, mixed-precision arithmetic, and enhanced data locality. By examining these methodologies, the discussion sheds light on how refinements in computational processes contribute to improved scalability and throughput, particularly in large-scale simulations. The complexities of heterogeneous memory systems are addressed in a dedicated section, which highlights the challenges of managing memory hierarchies in GPU-accelerated environments. Tools and frameworks like LLAMA are introduced as innovative solutions for optimizing memory layouts and minimizing data transfer overheads, offering practical insights into achieving efficient memory utilization. Finally, the report concludes with a summary of achievements and ongoing challenges, reflecting on the lessons learned and setting the stage for future advancements. It emphasizes the significance of these developments for the broader field of plasma simulations and high-performance computing, paving the way for exascale readiness and next-generation simulation capabilities. 2 Background on Accelerators and Heterogeneous Memory Using hardware accelerators in addition to standard central processing units (CPU) for either specific tasks, as a general purpose support or even as a replacement for most of the tasks usually performed by the CPU has been around for decades. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 7 Here, we specifically focus on accelerators used in scientific computing with a focus on accelerating numerical tasks. In recent years, the increased use of machine learning in connection with classical numerical tasks predominantly used for scientific simulations has broadened the spectrum of accelerator applications in scientific computing. Early applications of accelerators included the use of specific numerical co-processors that could execute certain numerical computations faster than standard CPUs [23]. Since the first decade of the century it became apparent that the progress in hardware used for graphics processing was fitting for scientific computing [12], starting a two decade long development in using general purpose graphics processing units (GPGPU ) for scientific computing now fueling the artificial intelligence revolution. Less prominent, yet in line with increasing parallelism in CPUs, cross-platform programming models [35, 41] and vectorization [31, 8] have seen increasing interest. Two of the major consequences of this development have been the creation of a variety of programming models, including vendor-based, closed source solutions [24] and open solutions [4, 13, 28] and the rise of deep memory hierarchies [22]. Both developments pose challenges to existing scientific simulation codes, as they usually have to use more than one programming model to run on a variety of state-of-the-art high-performance computer systems and are faced with the difficulties of complex data movement, hiding latencies and optimizing throughput. This has increased the need to optimize codes for specific accelerator hardware and to refactor codes for optimum data locality. As optimization is an additional burden to development teams, it has been observed that utilization of floating point performance in codes is often hampered by high latencies and low bandwidths in accelerator hardware, leading to many codes becoming memorybound [15, 16, 17]. The major part of this report will focus on GPGPU programming, especially on hardware from either Nvidia or AMD as the two major competitors in the field. When appropriate, we will refer to alternative subjects. With most of the European high performance computing systems in the Top 10 of supercomputers using accelerators from either Nvidia or AMD, developments for EPI accelerators is mainly deferred to WP2 and reported there. 3 Porting Status and Optimization Potential of the Codes This section examines the current state of the four codes, detailing their progress in utilizing GPU acceleration and heterogeneous memory. It reviews the specific optimization techniques applied to each code and identifies areas where further improvements are needed. 3.1 BIT Particle-in-Cell (PIC) Monte Carlo (MC) simulations are crucial for modeling plasmamaterial interactions, a critical component of fusion energy research due to the complex dynamics and the wide range of spatial and temporal scales involved. These simulations also incorporate atomic and collision processes. BIT1, a 1D3V PIC MC code, simulates D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 8 particle trajectories in a field governed by Maxwell’s and Poisson’s equations, enabling the study of plasma behavior with three-dimensional velocity components. Previous analyses of BIT1 identified performance bottlenecks, particularly in the particle mover function, which calculates particle trajectories and manages the movement of particles between grid cells and MPI ranks [39]. This stage is computationally intensive due to the need to process millions of particles [37]. Efforts have focused on enhancing performance by optimizing BIT1’s computational processes, particularly through asynchronous execution. BIT1 employs domain decomposition and MPI for parallelization, exchanging information at domain boundaries for tasks such as the smoother and Poisson solver. However, BIT1 currently does not utilize hybrid parallelization techniques, such as MPI+OpenMP or GPU offloading, relying solely on MPI. The code’s data layout stores particle positions and velocities within grid cells, transferring data between lists as particles move. Load imbalances arise in plasma simulations due to uneven particle concentrations across cells, reducing the efficiency of the particle mover function. Efforts to address these challenges aim to improve BIT1’s performance through asynchronous programming and optimization of its most demanding components. To identify optimal BIT1 porting opportunities, the focus has been on optimizing the particle mover function and analyzing its performance in large-scale plasma simulations. A test case was conducted to simulate neutral particle ionization in an unbounded, unmagnetized plasma, representative of conditions in fusion devices like ITER and DEMO. The scenario included electrons, D+ions, and Dneutrals, with ionization reducing neutral concentrations over time. The simulation employed a one-dimensional geometry with 100K cells, three plasma species, and 10M particles per species, totaling 30M particles across 200K time steps, excluding the field solver and smoother phases (shown in the diagram of [37]). The systems tested included Dardel, an HPE Cray EX supercomputer with AMD EPYC Zen2 processors and a 12 PB Lustre file system, running GCC 11.2.0 and Cray-MPICH; Vega, a EuroHPC petascale system featuring AMD EPYC 7H12 CPUs, Nvidia A100 GPUs, and Mellanox HDR100 InfiniBand, running GCC 12.3.0 and OpenMPI 4.1.2.1; and MareNostrum 5 (MN5), a EuroHPC pre-exascale system equipped with Intel Sapphire Rapids CPUs, Nvidia Hopper H100 GPUs, and ConnectX-7 NDR200 InfiniBand, operating with GCC 12.3.0 and OpenMPI 4.1.5. Figure 1: Hybrid BIT1 total simulation and optimized mover function using 16 and 64 ranks per node on Vega for 20000 times steps. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 9 3.1.1 Hybrid MPI and OpenMP/OpenACC BIT1 An in-depth investigation into BIT1’s performance was conducted, focusing on both intranode and inter-node testing to evaluate the impact of hybrid MPI and OpenMP/OpenACC approaches on execution time. The results, shown in Fig. 1, demonstrate that using MPI in combination with OpenMP or OpenACC significantly improves performance. With 16 MPI ranks, the hybrid MPI+OpenMP version reduced total simulation Figure 2: BIT1 optimized mover function - scaling up to 100 Nodes (12,800 MPI ranks) on Dardel for 200000 times steps. time by 53.4% (from 845.69 to 394.33 seconds) and mover function time by 69.3% (from 302.02 to 92.69 seconds). The MPI+OpenACC version further improved performance, reducing total simulation time by 57.0% to 355.85 seconds and mover function time by 70.3% to 89.61 seconds. Scaling to 64 MPI ranks, total simulation time dropped to 164.13 seconds and mover function time to 66.91 seconds, while adding OpenMP threads further reduced total simulation time to 42.07 seconds (a 74.4% reduction) and mover function time to 66.91 seconds (a 37.2% reduction). The best performance was achieved with the hybrid MPI+OpenACC version at 64 ranks, reducing both times to 32.32 seconds, representing an 80.3% reduction in total simulation time and a 51.7% reduction in mover function time compared to the baseline. Results from NJ and Vega show that BIT1 benefits from multicore CPUs through better core utilization via parallelization. Investigation into the scalability of hybrid BIT1 on CPUs, as shown in Fig.1, reveals a significant performance improvement as the number of MPI ranks increases from 2 (as reported in [37], D1.5 and D5.1) to 64, with noticeable improvements in both the total simulation and optimized mover function execution times, indicating that hybrid BIT1 scales effectively. To validate these findings, a second system, Dardel, was used for further investigation, specifically focusing on the optimized mover function with OpenMP (since OpenACC was not utilized). When scaling up to 100 nodes and 200,000 time steps on Dardel, our analysis demonstrates substantial performance gains with the hybrid MPI+OpenMP versions of BIT1 for both the mover function and total simulation. For the BIT1 mover function, as shown in Fig. 2, the MPI version with 1 node takes 260.85 seconds, while the hybrid MPI+OpenMP version improves this to 149.47 seconds, achieving a 42.7% speedup. Scaling to 10 nodes, the mover function time decreases to D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 16 OMP ACC ACC+ OMPX OMPX+ Timestep 4.7× 9.5× 14.8× 8.8× 13.3× Op Vlasov static 10.4× 162.7× 161.4× 226.9× 161.3× Op Vlasov dynamic 0.5× 5.8× 5.8× 5.5× 5.8× Op coll BGK 0.8× 38.3× 38.4× 32.0× 38.4× Ohm's law solver 2.3× 2.6× 3.4× 2.6× 3.2× Quasineutrality eq solver 0.9× 1.0× 1.5× 1.0× 1.5× Ampere's law solver 0.9× 1.0× 1.6× 1.0× 1.7× Op copy 1.4× 16.3× 16.6× 16.6× 16.5× Op axpy 1.0× 12.5× 12.6× 12.6× 12.6× Op linear combination 1.0× 14.8× 14.8× 14.9× 14.8× Speedups of each profiling region ACC ACC+ OMPX OMPX+ Timestep 2.0× 3.1× 1.9× 2.8× Op Vlasov static 15.6× 15.5× 21.7× 15.5× Op Vlasov dynamic 11.8× 11.8× 11.1× 11.8× Op coll BGK 50.7× 50.7× 42.3× 50.7× Ohm's law solver 1.2× 1.5× 1.2× 1.4× Quasineutrality eq solver 1.1× 1.6× 1.1× 1.6× Ampere's law solver 1.1× 1.8× 1.1× 1.8× Op copy 11.8× 12.0× 12.0× 11.9× Op axpy 12.5× 12.6× 12.6× 12.6× Op linear combination 14.6× 14.6× 14.8× 14.6× Speedups of each profiling region Figure 10: Speedups of each major profiling regions compared to the Fortran runtime (left) and the C++ runtime (right) on CPUs. The colormaps in the right figure are isolated in each rows. Scaling behaviour on the RZ dimension. The contrast difference between MCF reactors and experimental devices are in their poloidal size and the magnetic field strength. These two factors determine the numerical resolution on the poloidal cross section. In GENE-X, this translates to the number of points on RZ dimensions. Thus, scaling experiments on the RZ dimension is necessary for us to investigate the performance reliability of GENE-X in our target ITER case. Here we start the experiment from the same poloidal resolution as Table. 1 which has 212,387 RZ points the scaling up to a case with 2,621,778 RZ points which is around 12.3 times larger from where we start. In this experiment, we use the same number of nodes, GPUs, MPI processes as the previous case in Table. 1. The specifications of the other dimensions also stay constant and prescribed in Table. 2 The largest case here is estimated to use 45 GB of CUDA memory per GPU which is already 70.3% of the maximum capacity in MN5 ACC partition. This estimate only includes allocation on GENE-X side and does not consider allocation from inside the CUDA-powered multigrid GMRES solvers. The target ITER case is estimated to be around 20 times larger than this which means that we will run into GPU memory limit issue in the future. Table 2: Constant parameters for the scaling experiments on the RZ dimension with 20 nodes. Dimension Num. points MPI procs Num. points per procs φ16 4 4 v∥10 1 10 µ20 10 2 species 2 2 1 The scaling of the average elapsed time of both full Fortran run on CPU and OpenACC with CUDA-powered field solvers are shown in Fig. 11 by stacked area plot analogous D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 17 to the bar plots in Fig. 8. The stack of field solvers in red includes preprocessing time on GENE-X, i.e. matrix updates, while the dotted red stack shows the iterative solving time by the multigrid GMRES library PARALLAX and PAccX. The top figure of Fig. 11 shows near linear scaling behaviour of all profiling region in GENE-X. The preprocessing time of the field solvers is only a small fraction compared to the actual solving time. Similar to the result in the overview, Op Vlasov static is the bottleneck in the Fortran layer. This stencil operator is the largest and the most arithmetically complex in GENE-X. The bottom figure of Fig. 11 show super linear scaling behaviour of the total elapsed time. We can also see that the main contributor to such scaling behaviour is the field solvers. The performance of the CUDA-powered field solvers seems to be deteriorating as the problem size becomes larger in the RZ dimension. The multigrid GMRES solver only solves one poloidal plane at a time. This means that the kernel occupancy is solely determined by the number of RZ points. Despite of higher occupancy which commonly leads to better throughput, the CUDA kernel in the field solvers suffers in performance. This requires further investigation of the CUDA kernel in field solver with better profiler such as Nsight compute. 0.5 1.0 1.5 2.0 2.5 Number of points in RZ dimension ×106 0 500 1000 1500 2000 Avg. elapsed time of the time integration [s] Full Fortran Multigrid GMRES solvers Field solvers allreduce isend & irecv Misc. kernels Op coll BGK Op Vlasov dynamic Op Vlasov static 0.5 1.0 1.5 2.0 2.5 Number of points in RZ dimension ×106 0 1000 2000 3000 Avg. elapsed time of the time integration [s] OpenACC + CUDA Field solvers Multigrid GMRES solvers Field solvers allreduce isend & irecv Misc. kernels Op coll BGK Op Vlasov dynamic Op Vlasov static Figure 11: Average elapsed time of the time integration in GENE-X with respect of up-scaling the number of points in the RZ dimension. Full Fortran (top), combination of OpenACC and CUDA field solvers (bottom) are used. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 18 Figure 12: Weak scaling of PIConGPU on GH200 nodes taken from [11] for the JUPITER benchmark suite compared to other codes. Conclusion. The native compute kernels in GENE-X show remarkable performance improvements via OpenACC and OpenMP offload on GPUs. The results above are performed with CUDA-aware MPI in MN5. Optionally, we support NCCL for collective communication, i.e. allreduce, mainly for machine that is difficult to leverage CUDAaware MPI due to mixed-language mixed-vendor complication that we have in GENE-X. This complication occurs in Raven supercomputer in Max-Planck Computing and Data Facility. While we now fully support GENE-X time integration on GPU by combining directive-based approach, i.e. OpenACC and OpenMP offload, with CUDA and HIP for the field solvers, more performance improvements are still left to be desired in the field solvers. 3.3 PIConGPU PIConGPU [40] profits from cross-platform accelerator support using alpaka [41, 21, 20] for portable heterogeneous programming. alpaka uses a parallelisation hierarchy similar to that of CUDA and HIP, including grids, blocks, warps, threads and elements as levels in the parallel hierarchy, defining shared or non-shared memory access and parallelism for each level, see Fig. 14. As the levels are abstract they can be mapped to the actual parallelism hierarchy of a specific compute architecture depending on which programming model is selected as an alpaka backend [14]. In the course of the project SYCL was added to the backends, complementing backends such as CUDA, HIP/ROCm or OpenMP. In order to allow for extensive testing a complex CI/CD workflow to test backends for a variety of hardware platforms and compilers was set up as well as first kernels for continuous performance tests were developed. In addition to work on alpaka we also refactored and accelerated our parallel memory allocator / manager mallocMC built with alpaka to speed up small object allocation and manage device memory. For more details see section 5. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 19 101102103 # GCDs 109 1010 1011 1012 1013 FOM 90/10 LUMI strong scaling PIConGPU Figure 13: Strong scaling of PIConGPU on LUMI-G for the figure of merit as detailed in the PlasmaPEPSC project description. Green triangles indicate use of GPU-aware MPI while blue dots stand for standard MPI. We furthermore developed the LLAMA library for easy selection of memory layout for complex data types when switching compute architectures, see section 5. Using alpaka, PIConGPU is able to run the same simulation code on a large variety of devices, making adaption to specific EuroHPC systems a matter of recompilation. As an example, alpaka was used early in tests of GPU programming models for LUMI [20], showing superior performance compared to Kokkos and SYCL and on-par performance compared to native implementations. As a consequence, PIConGPU could present close to perfect weak scaling on both LUMI. Fig. 13 shows the strong scaling for PIConGPU’s figure of merit on LUMI-G, indicating that PIConGPU is on a good path towards reaching its goals. It was also selected for the JUPITER benchmark suite designed to estimate the performance of simulation codes for the upcoming EuroHPC exascale computer JUPITER as part of the JUPITER early access program JUREAP. Fig. 12 shows the weak scaling performance measured on JUWELS BOOSTER as part of the JUPITER early access program to estimate performance for the EuroHPC JUPITER system [11]. 3.4 Vlasiator Deliverable D5.1 described the approach taken in porting Vlasiator to GPU architectures, consisting of memory management, Vlasov solver porting, and field solver porting. The year of 2024 resulted in significant advancements in the operation of the Vlasov solver on GPU architectures. Memory management of the Vlasov solver and internal Vlasiator data structures was improved to the point where sample profiles were able to run several time steps without page faults of any data between the host and device. A significant performance bottleneck D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 20 was identified in use of many small kernel launches, resulting in device slowdowns. A major constituent of memory management, the block management of Vlasiator’s sparse velocity space, was completely re-designed to operate wholly on-device using parallel batch operations over the whole simulation domain. The majority of these improvements are detailed in a publication, submitted to the proceedings of the ASTRONUM-2024 conference, and available as a preprint [5]. The performance obtained in the code version described in this preprint is able to, for an artificial test case, exceed the performance of the CPU code path when comparing 1/4 CPU node against 1/4 CUDA GPU node on a HPC cluster. A detailed development plan exists for extending this performance to more generic use cases and improving performance overall. The porting of the field solver to GPU operations has proceeded as planned. Assemblylevel evaluation of loop operations accessing our field solver data motivated a more comprehensive overhaul of our field solver library FSgrid. A result is that FSgrid will be updated to provide a portable looping interface, resulting in portable code, which should provide optimized memory access patterns for both GPU and CPU use cases, without requiring a thorough re-write of the physics code describing the field solver equations. This is of especial importance when considering physical domain boundaries found in Vlasiator, such as the stepwise spherical inner boundary, requiring special management for different layers of the boundary domain. One major issue is an apparent performance gap between the CUDA and HIP versions of the code, running on NVIDIA and AMD hardware respectively. Implementing merged large overall kernels and pruning page faults have not yet succeeded in bridging this performance gap. This suggests that detailed analysis of AMD performance counters together is required in order to find the root cause. Current analysis indicates memory access patterns and cache utilization may play a role in this. 4 GPU Acceleration Strategies GPU accelerators, often called devices need a host, typically a CPU, to receive data for computation and initiate processes on the GPU. The process of copying data from host to device and back takes time and is limited by the latency and bandwidth between CPU and GPU. Such data transfer can be either evoked by the user or managed by a common memory interface which mirrors memory pages from the CPU to the GPU and vice versa. Copying data and starting processes from the host on the device is called offloading. Execution of processes on the device can be performed in parallel to processes on the CPU. At best, the waiting time for a full offloading cycle, meaning copying memory from host to device, executing code on the device and copying data back from device to host is at best spent doing computation in parallel on the CPU, e.g. for inter-node communication. In any case, the speed up of time to solution achieved by executing code on the accelerator should outperform execution on the CPU including the time needed for offloading. In the following we will discuss the many dimensions of using accelerator devices in PlasmaPEPSC and present lessons learned so far. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 21 4.1 Programming Interfaces Figure 14: Mapping of the alpaka parallel hierarchies Grid, Block, Warp, Thread and Element (right) to a multi-core CPU (left) and GPU (center). The two main vendors currently dominating the market for GPU accelerators are Nvidia, and to a smaller part, AMD. For two decades Nvidia has been providing the proprietary CUDA programming interface. AMD later introduced the HIP/ROCm interface that was heavily inspired by CUDA. In fact, AMD provides a HIP/CUDA translator, enabling the execution of many CUDA codes on AMD hardware without big changes. Non-vendor-specific programming models include OpenCL [9], which today is no longer competitive and has been surpassed by SYCL, directive-based models such as OpenMP [28] and OpenACC [27] or high-level APIs such as Intel’s OneAPI/SYCL [13]. While for many years only vendor compilers and tool chains were available to compile accelerator code, with the advent of [19] almost all major compilers, including Open Source compilers, understand accelerator code. Still, cross-platform code that can be executed efficiently on a variety of GPU platforms is sparse, as vendor-specific interfaces change rapidly and vendors invest in these rather than cross-platform solutions. Thus, cross-platform code usually adapts open standards (OpenMP, OpenACC, OpenCL). Maintaining these standards and providing hardware specific implementations requires significant resources and long-term support is not guaranteed. Thus, another approach is to develop light-weight programming interfaces that abstract across existing programming standards, providing a stable interface for portable parallel programming. Examples include Sandia’s Kokkos [35] and LLNL’s Raja [18] efforts. In PlasmaPEPSC the Helmholtz solution alpaka is used by PIConGPU, providing close-to-native performance for GPUs, CPUs and FPGAs. It provides mapping of a parallel hierarchy similar to CUDA including Grid, Block, Warp, Thread and Element to different accelerator back ends, see Fig. 14. 4.1.1 Multi-GPU Communication Interfaces On the lowest level, APIs such as CUDA provide multi-device pipelining in a single node. GPU-aware MPI [2, 25] now offers inter-node communication across nodes, multiple MPI D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 22 Infinity Fabric links: 4× AMD MI250X GPU1× AMD EPYC CPU NUMA 3 NUMA 3 NUMA 2 NUMA 2 partitioned L3 cache 205 GB/s 512GB -DDR4 GCD 7 GCD 6 GCD 0 GCD 1 36+36 GB/s GCD-CPU 50+50 GB/s GCD-GCD GCD0 GCD1 GCD6 GCD7 NUMA 1 NUMA 0 NUMA 2 NUMA 3 GCD 5 GCD 4 GCD 2 GCD 3 Figure 15: Architecture of a quad-GPU node of AMD MI250X, as featured in the LUMI supercomputer. processes sharing a single GPU and pipelining of MPI processes. Communication libraries like nccl [26] or rccl [3] complement MPI for multi-GPU, multi-node communication. Another library of interest is libfabric [1] which operates as a back end to e.g. MPI but also models for partitioned global address space approaches, which are not covered here. We provide an extensive characterization of the Infinity Fabric Interconnect [33]. This high-performance interconnect is present on multi-GPU AMD systems, and in particular is the intra-node interconnect used in the EuroHPC’s LUMI supercomputer, connecting together the four AMD MI250X GPUs on each node. Fig. 15 describes the complex mesh drawn by the Infinity Fabric Interconnect in such system. The first contributor to this complexity is the unique design of MI250X GPUs, which are composed, on a single board, of two independent Graphics Compute Dies (GCDs). A GCD is seen from the user perspective as a single GPU, with its own memory, and own compute units. In addition, the Infinity Fabric interconnect exposes four tiers of bandwidth, depending on the observed link. First, CPU-GPU Infinity Fabric links have a 36 GB/s bandwidth in each direction; in addition, GCD-GCD Infinity Fabric links have three tiers of bandwidth, respectively 200 GB/s per direction between same-GPU GCDs, 100 GB/s per direction for the pairs GCD2-GCD4 and GCD0-GCD6, and 50 GB/s per direction for all other pairs of connected GCDs. Furthermore, in this configuration, several hops might be required for data to transit within some pairs of GCDs, for example, between GCD5 and GCD6, one additional hop through GCD7 or GCD4 would be required. Those characteristics create a high-complexity mesh, which user might find trouble some to leverage directly. Thankfully, high-level interfaces, such as RCCL and GPU-aware MPI, abstract this complexity within a user-friendly set of API. From there, the question of performance and ability of those interfaces to leverage the full link bandwidth is raised. To answer this question, we evaluate the performance and crucial user-guided aspects of using those highlevel interfaces to leverage this interconnect. We further evaluate lower-level interfaces, such as explicit data movements with HIP, and GPU-side direct access to remote CPU memory. Explicit Movements. We first evaluate the CPU-GPU data movements, where we copy data from CPU memory to GPU memory using the dedicated hipMemcpy API. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 23 For this experiment, we evaluated two variants, both from pageable memory, simply allocated on the CPU-side with malloc, and pinned memory, allocated with hipHostMalloc. Our results show that movements from pageable memory achieved 71% of the theoretical bandwidth, while pinned memory achieves 79% of the theoretical bandwidth. This indicates quantitatively how memory pinning is essential to achieve the highest performance for data movements. For GPU-GPU peer-to-peer data movements, our experiments demonstrate the inability of the default configuration to leverage the full Infinity Fabric link bandwidth. Indeed, by default, explicit peer-to-peer data movements using hipMemcpyPeer rely on System Direct Memory Access (SDMA) Engines. This hardware feature is responsible for moving data in a way that copy can be overlap with GPU kernel execution. When using those SDMA engines, we observe a low utilization of the Infinity Fabric link, as low as 25% for data movements between two GCDs located on the same physical GPU. The use of SDMA engines can be disabled by setting the environment variable HSA ENABLE SDMA=0; however, setting this options removes the ability to overlap kernel execution with data movements. Therefore, a naturally-appearing alternative, which does not rely on SDMA engines, neither blocking on explicit data movements, is to use direct access from GPU kernel code. Direct access. We evaluate the performance of GPU-side kernel-level access to remote, CPU-located data. In this configuration, often referred to as Unified Memory, a single pointer is used to designate memory, which is located either on CPU or GPU memory. In our experiment, we allocate managed memory, using the hipMallocManaged API, and ensure an initial placement on CPU memory. We then use a GPU kernel, reading from this allocation. For this first experiment, we show that direct access over Infinity Fabric to CPU-located data achieves up to 70% of the theoretical Infinity Fabric link bandwidth, which is comparable to the level of performance obtained with explicit data movements. This shows that direct access from GPU kernel code to remote memory is a viable option for developers to use in their application in terms of performance, while removing the need for finely-tuned explicit data movements. We complete our analysis of GPU kernel-level access to CPU memory with an evaluation of the XNACK feature present in MI250X GPUs. This feature, used in combination with managed memory, migrates data to the GPU when accessed from kernel code, instead of performing a read over the Infinity Fabric Interconnect. This is aimed at improving data locality. Our experiments show that this approach exhibited as low as 8% of the theoretical CPUGPU Infinity Fabric link bandwidth. This is expected, as the overhead of migrating data on the first access can only be counterbalanced by repeated future access to this same data. As our benchmarking does not contain any data reuse, this benefit is not observed. However, XNACK feature can be potentially provide performance improvement when using unified memory in GPU-accelerated applications. In terms of peer-to-peer GPUGPU data movements, this options allowed reaching higher bandwidth than when using explicit data movements, generally on the level of 50% of the theoretical Infinity Fabric bandwidth. High-level Interfaces, RCCL and MPI. We evaluate GPU collective communications, where several GPUs execute collective operations. For this purpose, we compare the RCCL library, which is specialized for AMD GPUs, with GPU-aware MPI, which provides a generic API for all types of accelerators, and CPU-only systems. We use three types of collectives, one-to-all, where one GPU sends data to all others, all-to-one, D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 24 where all GPUs send data to a single one, and all-to-all, where all GPUs both send and receive data. Our results show that RCCL consistently outperforms GPU-aware MPI, with regards to latency for collective operations on 1 MB data size. We suggest that the multi-process approach used in MPI causes high overhead of inter-process communication, compared to the RCCL approach. Further study is however required for inter-node GPU collectives, as the behavior of the RCCL library and GPU-aware MPI might differ when compared to intra-node. PIConGPU makes use of NCCL, RCCL, GPU-aware MPI and libfabric either for inter-process communication or in-memory coupling of codes (see D4.2) on systems where these libraries are available. 4.2 Performance Optimisation Strategies Figure 16: Comparing the strategy to accelerate critical paths in a code (left) vs. using an abstraction layer to separate code implementation from accelerated code (right). Strategies for accelerating existing CPU code conventionally used critical path analysis to identify those parts in a code that dominated the overall execution time and could profit from the parallelization on the accelerator, see Fig. 16, left. However, while in theory this strategy sounds sensible, in reality one observes that instead optimizing the data flow of the application is key to achieve optimization. This means that offloading usually requires changing data formats and data flow paradigms of the code, resulting in major code changes all over the source code. Porting efforts then usually develop into porting more and more parts of the code to the accelerator, to minimize data movement. We have observed this when porting BIT to GPUs where time consuming parts of the code were ported one by one. A more time consuming strategy that usually results in better performance is to refactor the application code into user-facing code that provides a convenient interface to access all data and encapsulates all performance critical code via an API, see Fig. 16, right. This design is chosen, albeit in two distinct ways, by Vlasiator and PIConGPU. Another strategy for porting is to use high level libraries that provide optimizations for accelerators out of the box. This approach is heavily used by GENE/GENE-X. It has the advantage that the optimization is done by the library maintainers. Furthermore, D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 25 code can use the high level library instructions that typically abstract the underlying programming model and can provide portability across various platforms. A disadvantage is the strong dependency on these libraries, which must provide the performance needed by the application and be maintained for a long time. 4.2.1 GPU Matrix Processing Units Matrix processing units have been inspired by the rise of machine learning (ML) and artificial intelligence (AI) as yet another way of accelerating matrix-based operations. With growing interest in ML/AI the capabilities of these units have grown and added significantly to the overall compute power of GPUs. They are best for repetitive matrix operations or, more general, tensor operations, which are frequent in AI/ML. Similarly, some codes such as GENE/GENE-X use a lot of linear algebra operations from existing libraries, fast Fourier transforms and other matrix operations that benefit from vendorspecific acceleration through matrix processing units. Nvidia’s implementation of GPU matrix processing units, named Tensor Cores, have been widely studied in the literature. However, the study of AMD’s implementation of those matrix processing units, namely Matrix Cores, is still sparse in the literature. However, as EuroHPC’s LUMI supercomputer features such hardware unit, as part of the AMD MI250X GPU, evaluation of Matrix Cores is a crucial element for efficient codes. Therefore, we provide a first study of AMD Matrix Core [32], where we discuss the opportunities of programming those hardware units, from the lowest level of compilerlevel intrinsics, up to highest level of linear algebra libraries, such as GPU-accelerated BLAS. We associate this study of programming interfaces by a systematic performance evaluation, along with measurement of power efficiency. Low-level benchmarking. We first performed an evaluation of the low-level rocWMMA library, which directly maps C-level API calls to Matrix Cores instructions, with minimal to no overhead. This library only allows simple matrix operations, which are directly mappable to Matrix Cores abilities, therefore, the shape and datatypes of the matrices are strictly constrained. For this evaluation, we developed a micro-benchmarking approach, where we execute repeated Matrix Core instructions on all the available hardware units, for a significant amount of time. High-level programmability. The study of the highest-level rocBLAS library highlighted some important considerations with regards to mixed-precision GEMM operations. Using the singleand double-precision GEMM routines in rocBLAS, we were able to achieve a high percentage of floating-point throughput, relative to expected performance in the hypothesis of computations fully executed on Matrix Cores. Fig. 17 presents the results for various matrix sizes. We observe that the rocBLAS library achieves 90% and 100% of the throughput measured in our ideal microbenchmark setup, for DGEMM and SGEMM, respectively. This indicates that rocBLAS provides close-to-ideal utilization of Matrix Cores for singleand double-precision GEMM operations, without any low-level programming efforts. We demonstrate in our study that the performance gap between doubleand single-precision Matrix Core operations is induced by a high power consumption observed in double-precision, causing performance degradation to preserve hardware integrity. Therefore, it is advisable, both from the efficiency and raw performance perspectives, that single-precision is used when applicable, instead of double-precision. This D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 32 Figure 21: SNC-4 configuration of the Sapphire Rapids Figure 22: Latency on x86 Sapphire Rapids: DDR (left) and HBM(right) agement on the CPU. This approach would simplify handling both HBM and DDR while offering a unified framework for CPUs and accelerators. It is worth noting that SiPearl is also exploring additional frameworks to enhance data handling across memory tiers. Memory Exposure and Performance Considerations. The co-design of Rhea necessitates decisions about how memory is presented to users. The simplest approach is to expose one NUMA node for HBM and one for DDR. However, this abstraction obscures the hardware’s complexity, as multiple memory banks and controllers are aggregated into a single node. This simplification can negatively impact performance, particularly in scalar applications where the HBM targeted may be distant from the core in the Network-on-Chip (NoC). To maximize hardware efficiency, Sapphire Rapids (X86) architecture divides its chip into four quadrants, each paired with one HBM and one DDR. Although the design aggregates the two DDRs per quadrant into a single unit for simplicity (details Fig. 21), this configuration still provides insights into memory performance. For example, Fig. 22 illustrates memory access latency within and across quadrants for both HBM and DDR. Key observations from this architecture include: •HBM latency is approximately 30% higher than DDR latency. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 33 cores 5 10 20 40 DDR 409 sec 242 sec 136 sec 91 sec HBM 392 sec 201 sec 115 sec 72 sec speedup 4% 20% 19% 27% Table 4: Runtime impact of HBM/DDR for bit EBC.inp •Cross-quadrant memory access incurs a 15ns penalty for both HBM and DDR. •The worst-case penalty arises from accessing memory on the diagonal opposite quadrant, with a latency increase of up to 22%. These results underscore the importance of binding cores to the corresponding HBM/DDR pair within the same quadrant to achieve optimal performance. Since Rhea’s memory architecture is expected to be similar, this study provides valuable guidance for its co-design and memory exposure strategies. Early performance impact of HBM/DDR on PLASMA-PEPSC code Building on the considerations outlined earlier, initial experiments were conducted using two PLASMA-PEPSC codes: BIT1 and Vlasiator. For BIT1, it’s important to note that the extracted kernels are not ideal candidates for evaluating HBM and DDR performance. This is because these kernels are purely sequential, meaning they inherently favor DDR for two main reasons: (1) a single core cannot saturate the bandwidth of either HBM or DDR, and (2) DDR offers better latency. As a result, the experiments were carried out directly on the full BIT1 code. Consequently experiments have been performed directly on the BIT1 code. Table 4 summarizes the results of these experiments. The data shows that BIT1 achieves greater speedup as more cores are utilized: this highlights the benefits of using HBM for bandwidth bounded code. However, significant variability was observed during the tests, and further investigation is needed to understand the causes of this variability. For Vlasiator, Fig. 23 and Fig. 24 show the initial experiments conducted for both Magnetosphere 3D small and Ionosphere small test cases. In both cases, the configuration of 80 MPI / 1 OMP delivers the best speedup, achieving around 30%. However, switching to a configuration of 1 MPI / 80 OMP significantly improves overall performance, effectively eliminating the benefit of HBM. Additional experiments using more (hyper) threads followed this approach but failed to further enhance HBM’s performance gains. Detailed profiling revealed that the initial configuration uses twice as much memory as the alternative configuration. This increased memory usage transforms the test cases into bandwidth-bound scenarios. As a result, further investigations at larger scales are required. While evaluating these benchmarks, we also assessed the performance of the SYCL framework (see Delivery 2.4) on the HecBench benchmark. In this study, we compared the speedup achieved using HBM against DDR memory across 22 benchmarks that are compatible with all frameworks, ensuring consistency in comparison. All benchmarks D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 34 Figure 23: Impact of high bandwidth memory, presented as speedup for various MPI/OpenMP configurations on the Ionosphere small test case Figure 24: Impact of high bandwidth memory, presented as speedup for various MPI/OpenMP configurations on the Magnetosphere small test case D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 35 Framewok Minimum Median Mean Maximum OMPTarget 43% 101% 108% 303% ACPP omp,accelerated 83% 100% 116% 313% ACPP omp,library-only 92% 100% 108% 190% ACPP SSCP omp 92% 100% 117% 334% ACPP SSCP PoCL 92% 101% 116% 320% DPC++ PoCL 80% 101% 111% 247% DPC++ Native CPU 78% 99% 104% 169% Table 5: Speedup thanks to the HBM on various SYCL backend, targetting CPU, over a subset of 22 benchmarks coming from HecBenchmark were executed using 20 cores on Sapphire Rapids, which corresponds to the saturation point of the Stream bandwidth microbenchmark for this architecture. The results are summarized in Table 5. Several observations can be made from these results. First, certain benchmarks, notably the cross benchmark (a cross product of two 2D tensors), achieved an impressive 3x speedup. However, the AdaptiveCpp/LibOnly framework did not reach this level of performance. A closer examination of this framework revealed the use of nd range, which is strongly discouraged for this particular framework. At the opposite end of the spectrum, there was significant performance degradation observed for OMP-Target. Upon closer analysis, the affected benchmark was identified as atomicReduction, which repeatedly performs reductions in an atomic fashion. This highlighted the impact of latency, which is better managed by DDR in this case. Lastly, the overall performance gain was relatively modest on average. This is primarily because the HecBench suite is designed for general-purpose benchmarking, rather than being optimized for bandwidthor compute-bound kernels. Next Steps: This detailed analysis of Sapphire Rapids, along with the initial results, provides promising insights into the effects of NUMA exposure and MPI placement. SiPearl plans to refine these findings and develop more effective tools to optimize memory tiering and MPI placement, improving overall performance. 5.1.3 GPUDirect Storage Similar to remote direct data access that avoids routing memory copies between two GPUs through the host memory bus, GPUDirect Storage enables direct access to the file system without going through the CPU. In general, these operations are implemented by I/O libraries rather than directly by simulation code. One of the reasons for this is that the solutions are vendor-specific and thus not portable and optimization for a certain parallel file system can be cumbersome. The openPMD API for example mainly utilizes this via its back ends such as HDF5, as it does for GPUDirect RDMA for streaming in-memory workflows. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 36 5.2 Programming Interfaces for in-memory data Management in HMS In the following we focus on some solutions for in-memory data management of heterogeneous memory that work on different levels of the memory hierarchy. The list is neither complete nor exhaustive but covers various approaches towards managing heterogeneous memory management at varying complexity. The challenge is that vendor-specific solutions for memory management often do not perform well enough to be used for dynamical memory allocation during the simulation, especially for small object allocation. Moreover, memory handling can differ between vendors, thus abstractions that hide this complexity and can provide efficient allocation schemes are of big interest. We list a few approaches followed in PlasmaPEPSC below. 5.2.1 Concurrent memory allocation with mallocMC mallocMC is a policy-based C++ software framework developed by HZDR providing efficient, scalable and configurable on-device dynamic memory allocations orders of magnitude faster than what the vendor’s native implementations typically provide. The need for this kind of software is particularly high in PIC simulations where large numbers of particles are created, destroyed and moved around in every time step. Each particle being a small datum, the highly concurrent access to GPU memory being a shared resource brings challenges unseen on CPU. Consequently, standard CPU approaches to this problem cannot be ported and new algorithms have to be developed. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 37 100 101 102 All threads per block native malloc/free AMD Threadripper 1950X Nvidia A30 32 256 1k 2k 4k 16k 32k 256k Allocation size in bytes 100 101 One thread per block Speedup to native malloc/free Figure 25: mallocMC speed-up of one representative pattern of single-size allocations on an NVIDIA A30 as well as AMD Threadripper 1950X CPU compared to the vendor-native allocator. Larger is better. Please note the logarithmic axis. One sees that the AMD allocator performs better than the Nvidia one, while mallocMC provides either comparable or orders of magnitude better. mallocMC is maintained and developed as part of PLASMA-PEPSC in a completely open source fashion on Github.1In the latest iteration, it provides three different algorithms (plus the option for user-defined algorithms) for managing on-device memory within a kernel as well as configurable infrastructure around it. Preliminary benchmarks suggest a speed-up of up to 300x compared to native allocations on an NVIDIA A30 1https://github.com/alpaka-group/mallocMC D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 38 as can be seen in 25. A more thorough evaluation of its performance is under way. Its design makes it easy to maintain and extend such that keeping up with improvements on hardware and software side is fostered making it a future-proof and sustainable library to build on. At the time of writing, it has been powering PIConGPU for more than a decade already. While there is certainly some overlap of applicability, it should be stressed that mallocMC is on a significantly lower abstraction level than other solutions discussed here. For example, Umpire provides much higher-level memory management abstractions. mallocMC would rather correspond to a particular form of pool in Umpire. 5.2.2 Optimum memory layout and reduced representations with Llama Figure 26: LLAMA-provided mapping from user-defined data types (left) to an abstract index space (center) to a data layout in memory (right), from [10]. LLAMA, the library for Low-Level Abstraction of Memory Access, breaks with the C++ convention that the structure of a data type defines its layout in memory. While this paradigm had already been weakened by the introduction of the standard template library in C++, where data structure access and internal representation were partially disentangled, LLAMA uses standard C++ techniques to install a general, multi-dimensional index space that maps user-defined data types to this index space and maps this index space to a specific layout in memory, see Fig. 26. Examples for standard mappings already included in LLAMA are, array of struct, struct of array, array of struct of array. In addition, LLAMA provides user-defined mappings, working both at compile and run time. For mixed precision, mappings to reduced precision representations exist. Once data structures in code are annotated by LLAMA, switching to a given data layout is a matter of changing a single line of code. LLAMA was developed by HZDR together with CERN and Technical University Dresden and works independently of alpaka but is well suited for use together with alpaka, as alpaka maps processes to an abstract index space, making data parallel access equivalent to an index mapping operation. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 39 5.2.3 H2M H2M is a flexible, portable approach that works across different vendors and technologies to help manage data placement in systems with diverse types of memory. It allows programmers to specify requirements or provide hints—referred to as traits in the paper—about how data is accessed and used. The library uses these hints to efficiently allocate data across the available memory subsystems and can even move data dynamically at runtime. Additionally, they provide monitoring tools to make easier for programmers to choose the right traits for each piece of data by automatically generating memory access profiles and analyzing how data behaves, offering tailored recommendations for requirements or hints. While this library is suitable for early experiments, its flexibility and portable architecture come with the cost of storing each memory allocation on a single page. Consequently, this approach lacks of scalability and is underefficient in term of memory usage. Also this flexibility requires a lot of locking schemes on shared data-structures that alters performance. As a consequence, SiPearl is working on a more suited approach, combined with compiler directives that falls back into dedicated memory management library. Additionally, next step would be to have a look to UMF, a memory management framework provided by the oneAPI suite. 5.2.4 StarPU StarPU framework is a C++ wrapper that allows different ”codelets” to be implemented for heterogeneous architectures (EPI, RISC-V, FPGA, GPU) and run as tasks with datadependency awareness. Extracted kernels from work package 2 (WP2) can thus be implemented in an optimal way using dedicated compiler/language selection and specific optimisations. Calibration of codelets for performance models provide sizing of the buffers and timing under selected scheduling strategy. The main challenge in codelet design stays with data locality breakdown (cache coherency) for memory transfer between tasks and inter-node dependencies. Most of the StarPU development has been done with BIT1 code base for which memory reorganisation is underway. 5.2.5 Umpire Vlasiator and the Hashinator library [29, 30] are in the process of being ported to support memory allocation using the Umpire [6] portable allocation interface. Future work will allow a wider variety of memory allocators to be implemented. 6 Conclusion This report details strategies of porting Plasma-PEPSC codes to accelerator devices with a strong focus on GPUs. Various strategies, from accelerating critical paths in the codes to using high-level libraries such as GTENSOR or cross-platform parallel programming libraries such as alpaka are discussed. Particular emphasis is put on accelerator memory allocation, as vendor-specific solutions can be slow. Significant progress has been made across all Plasma-PEPSC codes in adapting to high-performance heterogeneous D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 40 computing systems. BIT1 porting to hybrid MPI and OpenMP/OpenACC frameworks has progressed significantly, achieving notable performance improvements and scalability on EuroHPC systems, with future efforts focusing on reducing GPU-CPU data transfer bottlenecks and exploring CUDA streams and hybrid approaches. GENE has reached production readiness with GTENSOR, while GENE-X has implemented GPU offloading for time integration, achieving substantial runtime reductions. Both codes continue to optimize AMD GPU performance, address field solver challenges, and prepare for largescale ITER simulations. PIConGPU has demonstrated strong scaling on platforms like LUMI and JUWELS BOOSTER through the alpaka library and enhanced memory management with mallocMC. Future efforts aim to explore using matrix processing units, and improve exascale performance with tests on JUPITER. Vlasiator has transitioned its Vlasov solver and memory structures to GPUs, enhancing stability and efficiency, though challenges remain in closing the performance gap between NVIDIA and AMD hardware. Work is on track to ensure the codes are well-prepared for the demands of exascale computing. We note that the project requires testing on systems that resemble the upcoming exascale systems to verify the results of research on the cutting-edge hardware. The project therefore would greatly benefit from the availability of the corresponding resources. D5.2: Algorithmic Changes and Memory Optimisation for Heterogeneous Systems 41 References [1] OpenFabrics Alliance. libfabric, 2024. https://ofiwg.github.io/libfabric/. [2] AMD. Gpu-aware mpi, 2024. https://gpuopen.com/learn/amd-lab-notes/ amd-lab-notes-gpu-aware-mpi-readme/. [3] AMD. Rccl, 2024. https://github.com/ROCm/rccl. [4] AMD. Rocm, 2024. https://rocm.docs.amd.com/en/latest/. [5] Markus Battarbee, Konstantinos Papadakis, Urs Ganse, Jaro Hokkanen, Leo Kotipalo, Yann Pfau-Kempf, Markku Alho, and Minna Palmroth. Porting the grid-based 3d+3v hybrid-vlasov kinetic plasma simulation vlasiator to heterogeneous gpu architectures, 2024. [6] David Beckingsale, Marty McFadden, Johann Dahm, Ramesh Pankajakshan, and Rich Hornung. Umpire: Application-focused management and coordination of complex hierarchical memory, May 2020. [7] Steven WD Chien and et al. sputnipic: an implicit particle-in-cell code for multi-gpu systems. In 2020 IEEE 32nd International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD), pages 149–156. IEEE, 2020. [8] Pawel Gepner. Using avx2 instruction set to increase performance of high performance computing code. Computing and Informatics, 36(5):1001–1018, Dec. 2017. [9] Khronos Group. Opencl, 2024. https://www.khronos.org/opencl/. [10] Bernhard Manfred Gruber, Guilherme Amadio, Jakob Blomer, Alexander Matthes, Ren´e Widera, and Michael Bussmann. Llama: The low-level abstraction for memory access. Software: Practice and Experience, 53(1):115–141, 2023. [11] Andreas Herten, Sebastian Achilles, Damian Alvarez, Jayesh Badwaik, Eric Behle, Mathis Bode, Thomas Breuer, Daniel Caviedes-Voulli`eme, Mehdi Cherti, Adel Dabah, Salem El Sayed, Wolfgang Frings, Ana Gonzalez-Nicolas, Eric B. Gregory, Kaveh Haghighi Mood, Thorsten Hater, Jenia Jitsev, Chelsea Maria John, Jan H. Meinke, Catrin I. Meyer, Pavel Mezentsev, Jan-Oliver Mirus, Stepan Nassyr, Carolin Penke, Manoel R¨ommer, Ujjwal Sinha, Benedikt von St. Vieth, Olaf Stein, Estela Suarez, Dennis Willsch, and Ilya Zhukov. Application-driven exascale: The jupiter benchmark suite, 2024. [12] Wolfgang H¨onig, Felix Schmitt, Ren´e Widera, Heiko Burau, Guido Juckeland, Mathias S. Muller, , and M. Bussmann. A generic approach for developing highly scalable particle-mesh codes for gpus. In Extended Abstract at 2010 Symposium on Application Accelerators in High Performance Computing (SAAHPC’10), Knoxville, USA, 2010, 2011. [13] Intel. Oneapi, 2024. https://www.intel.com/content/www/us/en/developer/ tools/oneapi/overview.html.