Plasma-PEPSC D4.2 Report on Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 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 D4.2 Report on Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms in Plasma-PEPSC Codes WP4: Extreme Data Analytics for Plasma Simulations Date of preparation (latest version): 31/12/2024 Copyright©2023 – 2027 The Plasma-PEPSC Consortium D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 2 DOCUMENT INFORMATION Deliverable Number D4.2 Deliverable Name Report on Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms in Plasma-PEPSC Codes Due Date 31/12/2024 Deliverable lead HZDR Authors Michael Bussmann (HZDR), Jeremy Williams (KTH), M˚ans Andersson (KTH), Carl-Martin Pfeiler (MPG), Franz Poeschel (HZDR), Diego Jimenez (MPG), Urs Ganse (UoH) and Yi Ju (MPG) Responsible Author Michael Bussmann (HZDR) E-mail: [email protected] Keywords Parallel I/O, In-situ, Compression WP/Task WP4/Task D4.2 Nature R Dissemination Level PU Final Version Date 31/12/2024 Reviewed by Leon Kos (UL), Stefan Costea (UL), Daniel Medeiros (KTH), Gabin Schieffer (KTH) and Luca Pennati (KTH) D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 3 DOCUMENT HISTORY Partner Date Comment Version KTH 03/03/2024 Skeleton version of the deliverable 0.1 MPCDF 20/11/2024 Dynamic in-situ technique in GENEX 0.2 KTH 26/11/2024 Compression method updates 0.3 KTH 28/11/2024 Parallel I/O and In-situ analysis updates in BIT1 0.4 UoH 29/11/2024 In-situ analysis techniques in Vlasiator 0.5 HZDR 02/12/2024 First draft 0.6 KTH 03/12/2024 Final draft updated for internal review 0.7 HZDR 12/12/2024 Revised draft after internal review 0.8 KTH 20/12/2024 Final cleanup for submission 1.0 D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 4 Executive Summary This report provides updates on the implementation of parallel I/O and its performance for selected Plasma-PEPSC codes. It dives deeply into several asynchronous and hybrid in-situ techniques for real-time data analysis and discusses lossless and lossy compression and reduced precision arithmetics for data storage, communication, and computation. These innovations aim to optimize data workflows, support large-scale simulations, and mitigate storage bottlenecks inherent to modern HPC systems. The report demonstrates substantial improvements in scalability, resilience, and computational throughput, essential for advancing plasma physics research in an exascale computing era. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 5 Contents 1 Introduction 6 2 Parallel I/O 7 2.1 New features of the openPMD-API . . . . . . . . . . . . . . . . . . . . . 8 2.2 BIT1 I/O Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 BIT1 I/O Costs Per Process . . . . . . . . . . . . . . . . . . . . . . . . . 11 3 In-situ data analysis 11 3.1 Improvements of in-situ data analysis in BIT1 . . . . . . . . . . . . . . . 12 3.2 Improvements of in-situ data analysis in GENE / GENE-X . . . . . . . . 12 3.2.1 Dynamic In-situ Techniques . . . . . . . . . . . . . . . . . . . . . 12 3.2.2 Implementation using MPI Sessions and DPP . . . . . . . . . . . 14 3.2.3 Evaluation of Dynamic In-situ Techniques . . . . . . . . . . . . . 16 3.3 Improvements of in-situ data analysis in Vlasiator . . . . . . . . . . . . . 18 3.4 Improvements of in-situ data analysis in PIConGPU . . . . . . . . . . . . 19 4 Compression Methods 22 4.1 Leveraging Compression Libraries for Communication . . . . . . . . . . . 22 4.2 Extending the availability of compression methods for scientific I/O . . . 25 4.3 Developing Statistical Methods for Plasma Simulation Data . . . . . . . 25 4.4 Reduced Precision Methods . . . . . . . . . . . . . . . . . . . . . . . . . 27 5 Conclusion 30 D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 6 1 Introduction This text reports on recent work in improving Parallel I/O, enabling in-situ data analytics workflows and data compression. High-Performance Computing (HPC) applications leveraging the continuously increasing peak performance of HPC systems can benefit data-intensive research. Plasma Physics simulations based on computationally expensive numerical methods profit immensely from this, producing more data than ever [33]. However, the storage capacity of HPC systems is restricted and this may limit the frequency of storing results for further analysis, which could limit scientific discovery. Another negative factor is that the input/output (IO) subsystem is also developing relatively slowly compared to the computational power and, thus, increasingly becoming a major bottleneck that limits an application’s performance. Conventionally, the results generated by an application are stored via the IO subsystem to the storage on HPC systems. This data is then read back by data analysis applications, again stressing the IO subsystems. Thus, enabling simulation codes to use the maximum bandwidth of parallel file systems is a first step in accelerating data-driven workflows. While compression can reduce the amount of data stored and potentially increase throughput, the time spent for compressing and decompressing data has to be taken into account as detailed in [23]. Compression can be used both in I/O and in-situ. In-situ techniques can help alleviate these problems by processing data directly at their place or on dedicated resources without writing and reading through the IO subsystem to and from the storage. Traditional in-situ techniques with static resource management can be categorized into three types depending on how the task interrupts the original application: in the synchronous approach (Fig. 1(a)), the in-situ task pauses the application until it has finished; in the asynchronous approach (Fig. 1(b)), extra resources execute the in-situ task, while the application continues its execution after having sent the data to these resources; in the hybrid approach (Fig. 1(c)), one part of the in-situ task pauses the original application, while separate resources execute the rest; it is a combination of the synchronous and asynchronous approach. In-situ workflows can avoid or reduce the IO throughput and storage, and proper choice among the three in-situ techniques can minimize the overhead brought by the in-situ techniques. Asynchronous and hybrid in-situ techniques are typically preferred for poorly scaling in-situ task [10, 24, 25]. However, they require separate resources. In many real-world scenarios, the in-situ processing will only start after a certain period. Therefore, these resources may stay idle for a significant amount of time. Strongly coupled in-situ techniques avoid data movement entirely, using in-place transformation of data [30, 32]. With dynamic resource management, one can allocate and release additional computing resources for in-situ tasks, avoiding leaving them idle, as shown in Fig. 1(d) and (e) compared to static resource allocations as shown in Fig. 1(b) and (c). In this document, we focus on the asynchronous and/or hybrid in-situ approaches, and thus refer to them when using the term in-situ, unless we explicitly specify the synchronous in-situ approach. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 7 Simulation In-situ task Data transfer (b) Application with static asynchronous in-situ task (a) Application with synchronous in-situ task (c) Application with static hybrid in-situ task Setup Main FinalIdle (d) Application with dynamic asynchronous in-situ task (e) Application with dynamic hybrid in-situ task Figure 1: Illustration of synchronous, asynchronous, and hybrid in-situ techniques with static resource management and asynchronous and hybrid in-situ techniques with dynamic resource management. Asynchronous or hybrid in-situ technique with dynamic resource management has less idle time than asynchronous or hybrid in-situ technique with static resource management. Strongly coupled in-situ techniques are not depicted, as they do not require data movement at all. The introduction of MPI Sessions [18] in the MPI Standard 4.0 [31] provides a new degree of flexibility and modularity in MPI, which could make MPI more suitable for elastic in-situ workloads [9]. A recent approach [12, 21, 20] introduced dynamic MPI Process interface extensions for MPI Sessions and an implementation based on Open MPI [13]. Complex asynchronous in-situ workflows, including reduction, compression and code coupling, are available via the openPMD library that accelerates I/O and provides inmemory coupling, as detailed in the following section. 2 Parallel I/O Increasing parallel I/O performance is crucial for •Writing checkpoints to increase resilience, •Writing high quality simulation data for later analysis •Reading simulation data from disk for further analysis or visualization. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 8 In the following sections we discuss recent I/O improvements to some Plasma-PEPSC codes. 2.1 New features of the openPMD-API The openPMD-api, a library utilized within the Plasma-PEPSC initiative by the simulation codes PIConGPU and BIT1, has recently been updated to version 0.16.0 [22]. This release introduces enhanced support for parallel I/O operations, further extending its compatibility with the ADIOS2 and HDF5 I/O frameworks. For HDF5, this includes support for the subfiling virtual file driver, a technique necessary for avoiding I/O bottlenecks in output procedures beyond a medium scale. More precise control over independent writing has been added in order to allow applications to use the most efficient I/O methods matching the parallel application logic. For ADIOS2, support for asynchronous background-writing has been added. The metadata aggregation logic has been greatly improved for extreme-scale applications with collective metadata logic (the usual case due to similar constraints imposed by HDF5). Support for the “joined array” dataset geometry has been added, reducing the number of collective operations needed for data preparation, and thus also reducing the API complexity for data preparation in simulation codes. Data files created with the openPMD-api via ADIOS2 can now optionally use an encoding that allows rapid inspection of data while it is being written, useful for inspecting partial data after a crash in parallel setups, enabling better diagnostics of errors. Support for ADIOS2 steps as an alternative to creating a new output file per simulation step has been improved, offering an I/O strategy with a far lower filesystem load for use cases with high-frequency output. An additional dataset geometry is available in a development version for using the metadata aggregation system in small-sized datasets, crucial for communication e.g. particle metadata in scalable streaming-I/O setups without requiring a fully connected communication mesh. For simplified debugging of parallel I/O, parallel JSON-based structure output has been added as an alternative to full output via HDF5/ADIOS2. 2.2 BIT1 I/O Improvements Introducing parallel I/O in Particle-in-Cell (PIC) Monte Carlo (MC) simulations is crucial, as it enables the efficient handling of data streams from multiple computational processes simultaneously. This parallel approach enhances throughput, reduces the time required for data storage and retrieval, and accelerates the pace of scientific discovery, thereby broadening the scope of simulations. Despite significant advances in computing capabilities, a persistent challenge remains in mitigating the performance bottleneck posed by I/O systems, which can hinder scientific progress and limit the potential of PIC MC simulations. BIT1, a 1D3V electrostatic PIC MC code, simulates plasma, impurity, and neutral transport in magnetic flux tubes of the plasma edge [39, 40, 41]. Despite its reduced dimensionality, BIT1 captures a broad range of kinetic processes, enabling novel plasma studies [38]. BIT1’s input is a compact file (1–3 kB), while the output workflows are defined by five critical parameters in Table 1. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 9 Table 1: Parameters for system diagnostics and simulation control. Parameter Description datfile Saves the system’s diagnostic snapshot at a specific time step. dmpstep Defines when the system’s current state, including particle data, is saved. mvflag Activates time-dependent diagnostics for plasma profiles and particle distributions. If >0, it sets the number of steps over which diagnostics are averaged. mvStep Sets the interval between time-dependent diagnostics. Controls how often the data is collected. last step Marks the final time step, saving the state to disk and ending the simulation. This ensures the completion of the simulation. As described in [45, 43], D4.1, and identified during development, BIT1’s original serial output supports up to 20,000 MPI processes effectively but struggles with largerscale simulations, where output becomes slow and prone to corruption. Addressing these limitations requires the implementation of parallel I/O methods to ensure scalability and data integrity for extensive HPC workloads. Figure 2: BIT1 Original File I/O Write Throughput, on Discoverer, Dardel and Vega CPU LFS, up to 200 Nodes, measured in GiB/s. [44] In [44], the openPMD-api parallel I/O library was integrated into BIT1 using the ADIOS2 backend BP4. The I/O performance of BIT1 was evaluated on the following three systems: •Discoverer: A petascale EuroHPC supercomputer featuring: –1128 compute nodes, D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 16 in-situ task(s), and the middle column shows the associated MPI process sets. Original Application Initialization: The original application processes are started with the mpirun command (2 processes in this example). After initializing an MPI Session, the processes create a common MPI communicator for the original application using the mpi://WORLD URI. Note that mpi://WORLD is a local alias for the URI PSet://1 containing the original application processes {p1, p2}. This is the MPI Session equivalent of using MPI COMM WORLD in the MPI World model. Subsequently, the original application uses this created communicator for the communication during the Original Application Steps. Resources Request: If additional resources are required, the resource manager is contacted via the extended MPI Sessions interface. To this end, an MPI PSet operation of type GROW is specified, where the URI PSet://1 is the input to the operation (black arrows in Fig. 8). The resource manager considers the specified process set operation in its scheduling algorithm and will eventually apply the advocated PSet operation by creating the two output PSets with URIs PSet://2 and PSet://3. These URIs are returned to the original application processes. PSet://2 refers to the newly created processes {p3, p4}(possibly on added resources) while PSet://3 refers to the union of PSet://1 and PSet://2. The newly created processes initialize an MPI Session and create a communicator for the in-situ task from the mpi://WORLD URI. Note, that here mpi://WORLD is a local alias for PSet://2 containing processes {p3, p4}. This is similar to disjoint MPI COMM WORLD communicators with MPI Comm spawn. Processes can determine if they were launched dynamically via the MPI info object associated with mpi://WORLD. To establish communication with the original processes a common communicator has to be created from PSet://3. To this end, the original processes publish the URI PSet://3 in the global dictionary associated with PSet://2, where the dynamically added processes can look up this information and participate in the communicator construction (purple arrows in Fig. 8). This common communicator can then be used for the Data Transfer from the original application to the in-situ task, e.g. after each application step. Application Finalization: To terminate application processes, these processes must disconnect from all communicators they are part of before finalizing the MPI Session. While MPI Comm disconnect is a collective and potentially synchronizing call on the processes within the communicator, MPI Session finalize is a local function. Therefore, the original application processes and the in-situ task processes can terminate independently. In this case, an MPI PSet operation of type SUB could be specified for a given PSet URI to release the resources associated with this PSet. However, in our implementation, we did not release resources prematurely, as the original application and in-situ tasks terminate almost simultaneously. 3.2.3 Evaluation of Dynamic In-situ Techniques All the experiments are performed on the Raven supercomputer at the Max Planck Computing and Data Facility (MPCDF) [2]. One CPU node contains two Intel Xeon IceLake-SP 8360Y processors, each with 36 cores and 256 GB RAM. We repeated the same experiments three times and only reported the average here because of negligible D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 17 Figure 9: Resource usage of GENE-X with in-situ image generation. The resource utilization for communication and in-situ tasks is barely visible and the dynamic asynchronous approaches (right two bars) have fewer idle resources differences among different runs. Dynamic Resource Management on HPC systems necessitates new optimization algorithms for global (dynamic) scheduling which is an active research area beyond the scope of our work. Since current production-level resource managers such as Slurm do not yet support dynamic scheduling we reserved enough resources for the original application and the in-situ task. Then, the application is started with less than the reserved resources, requesting in-situ resources dynamically over its runtime. We compute the resource usage cusing c= m X i ni×ti(1) with nithe occupied computing nodes, their respective occupancy times tiand mis the number of process groups. For our tests, we assume that resources are immediately available. While this might not always be the case on HPC systems with dynamic scheduling policies, our evaluation provides a baseline to understand the performance behavior of dynamic in-situ techniques. We executed a total of 200 simulation steps. From the 100th time step on, we performed one in-situ image generation every simulation step. We used four nodes for the D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 18 static synchronous and asynchronous in-situ approaches. We used two MPI processes on each node with 36 OpenMP threads per MPI process for the static synchronous approach. For the static asynchronous in-situ approach, we specified a CPU ID binding mask for each MPI process to correctly pin the threads to the process. We measured the shortest total execution time in our tests with the following configuration: two cores on each node are assigned to the two in-situ MPI processes, while the rest are to the simulation (two MPI processes and 35 cores per MPI process). Since image generation is not a computationally heavy task, we started the simulation with four nodes and added one additional node for the in-situ task in the basic dynamic approach. For the advanced dynamic approach, we started the simulation with four nodes, and the simulation automatically added one additional node for the in-situ task. As shown in Fig. 9, although the basic dynamic asynchronous approach avoids leaving the MPI process waiting for 100 simulation steps, it takes only a slightly shorter time than the static asynchronous approach, because the image generation is computationally cheaper than the simulation. The static asynchronous approach only uses two cores per node for image generation and avoids idle cores during image generation; the static synchronous approach requires no resources waiting for data, but it leaves the 35 cores per MPI processes idle when the image generation is executed; the default configuration of dynamic MPI processes added by MPI Sessions stays the same as the simulation, which leads to 35 cores per MPI process idle like the static synchronous approach. In this test case, the basic and advanced dynamic approaches use the same amount of resources. The total resource usage of the basic approach is slightly better due to a small difference in the communication time. 3.3 Improvements of in-situ data analysis in Vlasiator Vlasiator’s I/O subsystem is structured around the concept of a “data reducer”, with the idea of large-scale simulations as numerical experiments in mind. The most generic data reducers perform nothing but straightforward transformations of simulation data, such as writing electromagnetic field quantities or reduced velocity space moments into output files. With the increasing scale of global magnetospheric runs, it has become infeasible to write more than a handful of these quantities at full simulation cadence, and a more focussed in-situ analysis approach has been taken into action since the start of the Plasma-PEPSC project. Two examples of in-situ data analysis in Vlasiator that directly output science data are: •In [17], an in-situ data reducer running at high cadence during the simulation was used to obtain ionospheric ion precipitation spectra. By continuously classifying each phase space cell of the 6-dimensional magnetosphere simulation volume, it is possible to determine which parts of the near-Earth plasma population would follow the magnetic field lines into the ionosphere, where they can excite aurora and other ionospheric phenomena. The in-situ analysis process reduced the 6D simulation volume to a set of effectively 2D precipitation spectra, obviating the need to write significant amounts of phase space for post-processing. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 19 •In [36], massively parallel in-situ fieldline tracing with starting points in every simulation cell is employed to find the magnetic connectivity and topological structures in the magnetospheric magnetic fields. In this particular publication, the mechanism was used to find flux ropes – wound structures of magnetic field lines. The same mechanism can also be employed to determine which parts of the simulation is magnetically connected to Earth, meaning that plasma phenomena in those regions can directly impinge onto Earth’s atmosphere. Performing this step in-situ, using full electromagnetic field data and high-order derivative information, is both more precise and faster than a postprocessing approach. 3.4 Improvements of in-situ data analysis in PIConGPU The year 2024 has seen an exploration into practical solutions for performing in-situ application coupling for analysis on the raw data produced by PIConGPU at the full scale of Oak Ridge Leadership Computing Facility (OLCF) Frontier. In-situ hereby refers to a data strategy where data is exchanged from one application’s memory to another application’s memory, possibly located on another compute unit, without using permanent storage for the task. This strategy is chosen for a multitude of reasons: 1. Integrating all parts of a scientific data processing pipeline into a single application is often not possible due to vastly different software environments. The example at hand is a data processing pipeline consisting of the accelerated C++ simulation code PIConGPU and a PyTorch-based machine learning model for inferring the relation between the phase space characteristics and the particle trajectories causing the phase space. 2. Integrating all parts of a scientific data processing pipeline into a single application is often not desirable due to vastly different scaling properties of heterogeneous stages in a data processing pipeline. Here, PIConGPU scales up to full system size while the used PyTorch code scales up to only a medium portion of the system before being overwhelmed by collective communication operations. 3. Using the parallel filesystem is generally no longer possible on recent Exascale compute systems due to limited I/O bandwidth and capacity in relation to the amounts of raw data produced by accelerated simulation codes. 4. An alternative should allow flexible communication patterns without locking scientists into a tightly specified application setup, as techniques like shared memory or node-local SSDs for data exchange inevitably would. Communication must hence be in-transit, making use of HPC network interconnects. PIConGPU was tested against a synthetic data consumer to test scalability of Streaming I/O via the OLCF Frontier’s Slingshot network at medium to full scale, using different implementations to access the Slingshot network for remote direct memory access (RDMA). The measured performance results are seen in Fig. 10. The PIConGPU Kelvin–Helmholtz instability (KHI) simulation is used as a reference, producing 5.86 GB of particle data per compute node and time step. At each scaling run D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 20 Figure 10: Parallel total throughput for streaming at full scale on Frontier using a synthetic benchmark built on PIConGPU KHI. An obvious outlier result was removed for the libfabric backend at 8192 nodes. (a) Using the libfabric data plane and the CXI provider to access the Slingshot network at a lower level. (b) Using the MPI data plane based on MPI Open port(). The perceived parallel throughput reaches around 20TB/s for libfabric, and from 20 - 30 TB/s for MPI. (weak scaling), five time steps are sent from PIConGPU to the no-op consumer, which then measures the time needed for loading the data. The perceived parallel throughput is computed based on this measured time and the global data size, i. e. the previously mentioned 5.86 GB per node multiplied with the number of nodes. While including communication overhead, this perceived parallel throughput has been shown in [10] to be a close approximation for the real throughput in this kind of setup. The benchmarks depicted in Fig. 10 demonstrate the feasibility of this I/O strategy for full-scale workflows. Depicted are boxplots of all single measurements. Most notably, we observe a maximum perceived parallel throughput of 20 - 30 TB/s which compares outstandingly against the 10 TB/s bandwidth of the parallel Orion filesystem, a scaling limit which we circumvent and exceed by not using the filesystem for intermediate data. Furthermore, these results can also compare well against the 35 TB/s aggregate write bandwidth of the SSDs installed locally in the compute nodes [3]. For the libfabric-based benchmarks, we initially used an implementation that enqueued all read operations to the network at once and waited for replies (labeled 4096∗ in Fig. 10). While these results signify the best achieved per-node throughput of 3.5 - 4.7 GB/s at 4096 compute nodes (a total of 14.1 - 18.8 TB/s), it turned out that this strategy did not scale to the full system, so we employed an alternative that enqueued operations in batches of 10 operations. While this turned out to scale to the full system, it came at a notable cost in performance, visible in a per-node throughput of 1.9 - 2.6 GB/s at 9126 compute nodes (a total of 16.5 - 23.0 TB/s). Conversely, the MPI data plane yields a per-node throughput from 2.6 - 3.7 GB/s at 4096 compute nodes (a total of 10.5 - 14.9 TB/s) to 2.4 - 3.3 GB/s at 9126 nodes (a total of 21.4 - 29.5 TB/s),the best achieved system-wide parallel throughput from this experiment. We conclude that the MPI data plane brings default good performance, while the libfabric data plane’s lower-level control can bring performance improvements but at the cost of necessary fine-tuning beyond this study’s scope. Since only a single instance of D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 21 the reader code was running per node, its throughput (1.9 - 4.7 GB/s, across all cases) compares against the max possible throughput of a single HPE Slingshot NIC at 25 GB/s. This implies that further speedup can be achieved by parallelizing the reader in the actual PIConGPU+MLapp pipeline. The benchmarks suggest that the MPI data plane be used, which can be explained by leveraging system-specific performance tuning in the MPI implementation. All regular measurements range between 1.2 s - 3.2 s. Future investigations will include data aggregation patterns for streaming I/O in order to support the different sizes at which both ends of the data processing pipeline scale, as well as increased per-node parallelism for making better use of the network access reachable from processes executing on a compute node. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 22 4 Compression Methods The four flagship codes enable simulations that require a large amount of data to be handled during and after simulations. In this section, we describe how lossless and lossy compression can be used to speed up simulations, increase througput and reduce storage usage after simulations. During simulations, a limiting factor for performance is communication and data movement. Fast intra-node memory connectivity and network-based inter-node communication is crucial for the scalability of the code. The compression of data that is communicated between parallel instances of a code can improve commutation and the time to solution. A complete analysis must include the time for compression and decompression to estimate overall speed up due to compression. One of the main questions for such methods is to determine when communication with compression and decompression is faster than communication with non-compressed data. Multiple compression libraries have been evaluated and ZFP was chosen for the deeper evaluation, which is presented in section 4.1. During and after a successful simulation, a large amount of storage is required to store the particle and field information. This is a growing problem in the field as larger simulations with better resolution and more correct representation of physical phenomena generate more data. Storage is both costly and energy-intensive, presenting scientists with a challenging dilemma: either retain large volumes of data at significant expense or risk deleting data that is essential for preservation and accessibility within the scientific process. To enable the use of compression libraries in I/O the effort to integrate it in openPMD-api is described in section 4.2. These large simulations also require resilience measures, as interruptions that lead to canceled simulations are expensive (scientist’s time as well as core hours). Therefore, in I/O tasks such as check-pointing and long-term storage lossy compression is advantageous. This task has explored Gaussian Mixture Models (GMM) as a Machine Learning approach for the compression of particle data described in section 4.3 mainly intended to be used in check-pointing but can also be useful for post-processing and analysis of large simulations. Compression can be used to lower the memory consumption of a simulation, improving memory-bound algorithms. One way of compressing data is to represent it with a lower accuracy floating-point representation. This work focuses on experimental validation of analytic error estimates for spatial discretization and its required floating-point precision in finite difference calculations. This has been performed with a library designed to do floating-point calculations with arbitrary precision described in section 4.4. The study has concluded that some simulations waste memory and compute by over-resolving lowresolution problems. 4.1 Leveraging Compression Libraries for Communication As mentioned on the previous deliverable: D4.1 Initial Report on Improving Parallel I/O, Enabling In-situ and Compression, Max Planck Gesellschaft (MPG), through Max Planck Computing and Data Facility (MPCDF) and Max-Planck Institut f¨ur Plasmaphysik (MPG IPP), is involved in a German-funded project titled Data Reduction for Exascale Application in Fusion Research (DaREXA-F). Within this project, MPCDF D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 23 Figure 11: Measuring ZFP compression/decompression overhead with respect to uncompressed communications. Given a simple send-receive communication scheme within pairs of ranks holding their data in GPU device buffers, compressed communication is viable for data exchanges larger than 107bytes. has been exploring lossy compression methods with a focus on inter-node communication as part of distributed computations. In particular, after evaluating lossy compression libraries like ZFP [29, 7], MGARD [4, 27, 16] and SZ [6, 37, 26, 28] for performance and data compressibility, ZFP was chosen as the most performant and easy to adopt given GENE’s complex numerical data, comprising both real and imaginary components. One of the main challenges in using in-simulation data compression for communication is ensuring the compression and decompression processes are fast enough to minimize overhead, ideally making it negligible. Fig. 11 displays an evaluation performed with the ZFP compressor where GPU-based data is compressed before it is communicated and decompressed on the receiving end. Given different values of Compression Factor (CF) = Uncompressed Size Compressed Size , we can see that costs of compression overhead make it viable only to use compression for messages larger than approximately 10 to 12 MB. However, as shown in Table 2, for a target scenario of an ITER-like problem size and given a reasonable parallelization in GENE (assuming roughly 8640 GPUs are used), the biggest message exchanges would be in the order of 7.5 MB. We have developed a scheme that enables us to hide some of the compression and decompression costs by overlapping these operations with the communication of compressed data to reduce the message size threshold at which compressed communications provide performance returns. To do this, a modified version of the ZFP library that supports encoding and decoding kernels to be executed on user-defined GPU streams is used in hand with the NVIDIA Collective Communications Library (NCCL), which provides pointto-point operations for data exchange across devices streams. This overlapped scheme has been introduced into GENE’s v-direction boundary exchange operation in which the D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 24 Table 2: Expected dimensions of a message size in the exchange v 6d operation in global-mode GENE for an ITER-like scenario of burning-plasma. GENE buffers in this exchange operation are double precision complex numbers. Dimensions Size (Minimal ITER-like sizes) MPI Ranks sbuf shape Message Size nx0 1024 4 256 nky0 96 3 32 nz0 24 12 6 ∼7.5 MB nv0 32 4 2 nw0 25 5 5 n spec 3 3 1 data is then compressed and decompressed during each neighbor exchange. We have evaluated this new communication scheme using a Cyclone-Base case in GENE and ran experiments on the GPU partition of the Raven supercomputer at the MPCDF [2]. Each Raven GPU-accelerated node is comprised of 4 NVIDIA A100 GPUs, interlinked with NVLink 3 and the network connecting GPU nodes is a Mellanox Infiniband HDR200. Table 3: Dimensions of a message size in the exchange v 6d operation in the Cyclone-Base Case. Dimensions Size (Minimal ITER-like sizes) MPI Ranks sbuf shape Message Size nx0 96 1 96 nky0 28 1 28 nz0 32 4 12 ∼7.8 MB nv0 32 8 2 nw0 8 1 8 n spec 1 1 1 Table 3 shows the test case dimensions and the chosen parallelization for evaluation. The message size is close to the expected target ITER scenario and the parallelization is done in such a way that allows us to evaluate the proposed scheme under network saturation conditions, where the benefits of data compression would be useful for exascale simulations. Fig. 12 shows the mean time per timestep for the test case described in Table 3. Normal execution refers to a GENE simulation where no compression is used. This is also the baseline reference for performance. Given different levels of network usage in Raven, normal execution shows more performance variability than the compressed experiments. We use speedup ranges to quantify the effects of compression over the execution time of each time step in GENE and observe up to a 16% performance improvement when communication is done with a compression factor of 8. Furthermore, main quantities of interest resulting from compressed runs of the simulation like the spatially averaged estimated means for electrostatic heatflux and parallel velocity space moment are within [0.31% −10.80%] of the expected value for the former and [1.51% −2.09%] for the latter, highlighting the potential for lossy compression in GENE. In the context of this project, we will evaluate this proposed overlapped scheme also in other applications, in particular D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 25 Figure 12: Speedups resulting from proposed overlapped compression scheme with ZFP Vlasiator, where previous studies [11] have shown that this is a communication-bound code where message sizes are relatively large and involve high dimensionality data. 4.2 Extending the availability of compression methods for scientific I/O The openPMD-api, used currently among the simulation codes pertaining to PlasmaPEPSC by PIConGPU and BIT1, has so far supported compression options only via the operator mechanic of ADIOS2, but not yet via the filter mechanic of HDF5. While both these techniques are not restricted to compression methods, they offer a wide range of established compression operators, open to third-party contributors. Options include e.g. GZip, BZip2, Blosc2, SZ or ZFP. The recent release 0.16.0 of the openPMD-api [22] adds targeted configuration for HDF5 chunking, a necessary prerequisite for compression. A usable prototype for compression via HDF5 in the openPMD-api is now available as a development version1, thus adding to the I/O flexibility of simulation codes building on this I/O solution. 4.3 Developing Statistical Methods for Plasma Simulation Data There are multiple strategies for lossy compression based on Machine Learning and statistical methods, some described in deliverable 4.1. The main idea for using GMM in the compression of PIC data comes from the insight that in most of the simulations the velocity distribution function does not greatly deviate from a Maxwellian or a Gaussian. Therefore we can approximate the particle velocities as a Velocity Density Function (VDF) modeled as a linear combination of Gaussian distributions, which in a compressed state can be reduced to a set of mean values and variance/co-variance matrices (in higher 1https://github.com/openPMD/openPMD-api/pull/1644 D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 32 [13] Edgar Gabriel, Graham E. Fagg, George Bosilca, Thara Angskun, Jack J. Dongarra, Jeffrey M. Squyres, Vishal Sahay, Prabhanjan Kambadur, Brian Barrett, Andrew Lumsdaine, Ralph H. Castain, David J. Daniel, Richard L. Graham, and Timothy S. Woodall. Open mpi: Goals, concept, and design of a next generation mpi implementation. In Dieter Kranzlm¨uller, P´eter Kacsuk, and Jack Dongarra, editors, Recent Advances in Parallel Virtual Machine and Message Passing Interface, pages 97–104, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg. [14] K. Germaschewski and B. Allen. GTensor is a multi-dimensional array C++14 header-only library for hybrid GPU development. [15] K. Germaschewski, B. Allen, T. Dannert, M. Hrywniak, J. Donaghy, G. Merlo, S. Ethier, E. D’Azevedo, F. Jenko, and A. Bhattacharjee. Toward exascale wholedevice modeling of fusion devices: Porting the GENE gyrokinetic microturbulence code to GPU. Physics of Plasmas, 28(6):062501, 06 2021. [16] Qian Gong, Jieyang Chen, Ben Whitney, Xin Liang, Viktor Reshniak, Tania Banerjee, Jaemoon Lee, Anand Rangarajan, Lipeng Wan, Nicolas Vidal, et al. Mgard: A multigrid framework for high-performance, error-controlled data compression and refactoring. SoftwareX, 24:101590, 2023. [17] Maxime Grandin, Thijs Luttikhuis, Markus Battarbee, Giulia Cozzani, Hongyang Zhou, Lucile Turc, Yann Pfau-Kempf, Harriet George, Konstantinos Horaites, Evgeny Gordeev, Urs Ganse, Konstantinos Papadakis, Markku Alho, Fasil Tesema, Jonas Suni, Maxime Dubart, Vertti Tarvus, and Minna Palmroth. First 3d hybridvlasov global simulation of auroral proton precipitation and comparison with satellite observations. Journal of Space Weather and Space Climate, 13:20, 2023. [18] Daniel Holmes, Kathryn Mohror, Ryan Grant, Anthony Skjellum, Martin Schulz, Wesley Bland, and Jeffrey Squyres. Mpi sessions: Leveraging runtime infrastructure to increase scalability of applications at exascale. pages 121–129, 09 2016. [19] Dominik Huber. Prototype, 2023. [20] Dominik Huber, Martin Schreiber, Martin Schulz, Howard Pritchard, and Daniel Holmes. Design principles of dynamic resource management for high-performance parallel programming models, 2024. [21] Dominik Huber, Maximilian Streubel, Isa´ıas Compr´es, Martin Schulz, Martin Schreiber, and Howard Pritchard. Towards dynamic resource management with mpi sessions and pmix. In Proceedings of the 29th European MPI Users’ Group Meeting, EuroMPI/USA’22, page 57–67, New York, NY, USA, 2022. ACM. [22] Axel Huebl, Franz Poeschel, Fabian Koller, Junmin Gu, Michael Bussmann, JeanLuc Vay, and Kesheng (John) Wu. openPMD-api: C++ & Python API for Scientific I/O with openPMD, June 2018. [23] Axel Huebl, Ren´e Widera, Felix Schmitt, Alexander Matthes, Norbert Podhorszki, Jong Youl Choi, Scott Klasky, and Michael Bussmann. On the scalability of data D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 33 reduction techniques in current and upcoming hpc systems from an application perspective. In Rio Yokota, Julian M. Kunkel, Michela Taufer, and John Shalf, editors, High Performance Computing - ISC High Performance 2017 International Workshops, DRBSD, ExaComm, HCPM, HPC-IODC, IWOPH, IXPUG, P3MA, VHPC, Visualization at Scale, WOPSSS, Revised Selected Papers, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), pages 15–29. Springer Verlag, 2017. Publisher Copyright: ©Springer International Publishing AG 2017.; 32nd International Conference on High Performance Computing, ISC High Performance 2017 ; Conference date: 1806-2017 Through 22-06-2017. [24] Yi Ju, Mingshuai Li, Adalberto Perez, Laura Bellentani, Niclas Jansson, Stefano Markidis, Philipp Schlatter, and Erwin Laure. In-situ techniques on gpu-accelerated data-intensive applications. In 2023 IEEE 19th International Conference on eScience (e-Science), pages 1–10. IEEE, 2023. [25] Yi Ju, Adalberto Perez, Stefano Markidis, Philipp Schlatter, and Erwin Laure. Understanding the impact of synchronous, asynchronous, and hybrid in-situ techniques in computational fluid dynamics applications. In 2022 IEEE 18th International Conference on e-Science (e-Science), pages 295–305. IEEE, 2022. [26] Xin Liang, Sheng Di, Dingwen Tao, Sihuan Li, Shaomeng Li, Hanqi Guo, Zizhong Chen, and Franck Cappello. Error-controlled lossy compression optimized for high compression ratios of scientific datasets. In 2018 IEEE International Conference on Big Data (Big Data), pages 438–447. IEEE, 2018. [27] Xin Liang, Ben Whitney, Jieyang Chen, Lipeng Wan, Qing Liu, Dingwen Tao, James Kress, David Pugmire, Matthew Wolf, Norbert Podhorszki, et al. Mgard+: Optimizing multilevel methods for error-bounded scientific data reduction. IEEE Transactions on Computers, 71(7):1522–1536, 2021. [28] Xin Liang, Kai Zhao, Sheng Di, Sihuan Li, Robert Underwood, Ali M Gok, Jiannan Tian, Junjing Deng, Jon C Calhoun, Dingwen Tao, et al. Sz3: A modular framework for composing prediction-based error-bounded lossy compressors. IEEE Transactions on Big Data, 9(2):485–498, 2022. [29] Peter Lindstrom. Fixed-rate compressed floating-point arrays. IEEE transactions on visualization and computer graphics, 20(12):2674–2683, 2014. [30] Alexander Matthes, Axel Huebl, Ren´e Widera, Sebastian Grottel, Stefan Gumhold, and Michael Bussmann. In situ, steerable, hardware-independent and data-structure agnostic visualization with isaac. Supercomputing Frontiers and Innovations, 3(4):30–48, Oct. 2016. [31] Message Passing Interface Forum. MPI: A Message-Passing Interface Standard Version 4.0, June 2021. [32] Felix Meyer, Benjamin Hernandez, Richard Pausch, Ren´e Widera, David Groß, Sergei Bastrakov, Axel Huebl, Guido Juckeland, Jeffrey Kelling, Matt Leinhauser, D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 34 David Rogers, Ulrich Schramm, Klaus Steiniger, Stefan Gumhold, Jeff Young, Michael Bussmann, Sunita Chandrasekaran, and Alexander Debus. Hardwareagnostic interactive exascale in situ visualization of particle-in-cell simulations. In Proceedings of the Platform for Advanced Scientific Computing Conference, PASC ’23, New York, NY, USA, 2023. Association for Computing Machinery. [33] Dominik Michels, Andreas Stegmeir, Philipp Ulbl, Denis Jarema, and Frank Jenko. Gene-x: A full-f gyrokinetic turbulence code based on the flux-coordinate independent approach. Computer Physics Communications, 264:107986, 2021. [34] Dominik Michels, Philipp Ulbl, Wladimir Zholobenko, Thomas Body, Andreas Stegmeir, Thomas Eich, Michael Griener, Garrard D Conway, Frank Jenko, ASDEX Upgrade Team, et al. Full-f electromagnetic gyrokinetic turbulence simulations of the edge and scrape-off layer of asdex upgrade with gene-x. Physics of Plasmas, 29(3), 2022. [35] Truong Nguyen, Guangye Chen, and Luis Chac´on. An adaptive EM accelerator for unsupervised learning of Gaussian mixture models. arXiv preprint arXiv:2009.12703, 2020. [36] Yann Pfau Kempf, Konstantinos Papadakis, Markku Alho, Markus Battarbee, Cozzani Giulia, Lauri P¨ank¨al¨ainen, Urs Ganse, Fasil Kebede, Jonas Suni, Konstantinos Horaites, Maxime Grandin, and Minna Palmroth. Global evolution of flux transfer events along the magnetopause from the dayside to the far tail. Annales Geophysicae, 2024. [37] Dingwen Tao, Sheng Di, Zizhong Chen, and Franck Cappello. Significantly improving lossy compression for scientific data sets based on multidimensional prediction and error-controlled quantization. In 2017 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 1129–1139. IEEE, 2017. [38] D Tskhakaya. One-dimensional plasma sheath model in front of the divertor plates. Plasma Physics and Controlled Fusion, 59(11):114001, 2017. [39] D Tskhakaya and S Kuhn. The magnetised plasma-wall transition: Theory and pic simulation. Contributions to Plasma Physics, 44(5-6):564–570, 2004. [40] David Tskhakaya and et al. Optimization of pic codes by improved memory management. Journal of Computational Physics, 225(1):829–839, 2007. [41] David Tskhakaya and et al. Pic/mc code bit1 for plasma simulations on hpc. In 2010 18th Euromicro Conference on Parallel, Distributed and Network-based Processing, pages 476–481. IEEE, 2010. [42] Jeremy J Williams, Stefan Costea, Allen D Malony, David Tskhakaya, Leon Kos, Ales Podolnik, Jakub Hromadka, Kevin Huck, Erwin Laure, and Stefano Markidis. Understanding the impact of openpmd on bit1, a particle-in-cell monte carlo code, through instrumentation, monitoring, and in-situ analysis. In European Conference on Parallel Processing. Springer, 2024. D4.2: Improving Parallel I/O, and Enabling In-Situ and Compression Algorithms 35 [43] Jeremy J Williams, Felix Liu, David Tskhakaya, Stefan Costea, Ales Podolnik, and Stefano Markidis. Optimizing bit1, a particle-in-cell monte carlo code, with openmp/openacc and gpu acceleration. In International Conference on Computational Science, pages 316–330. Springer, 2024. [44] Jeremy J Williams, Daniel Medeiros, Stefan Costea, David Tskhakaya, Franz Poeschel, Ren´e Widera, Axel Huebl, Scott Klasky, Norbert Podhorszki, Leon Kos, et al. Enabling high-throughput parallel i/o in particle-in-cell monte carlo simulations with openpmd and darshan i/o monitoring. In 2024 IEEE International Conference on Cluster Computing Workshops (CLUSTER Workshops), pages 86– 95. IEEE, 2024. [45] Jeremy J Williams, David Tskhakaya, Stefan Costea, Ivy B Peng, Marta GarciaGasulla, and Stefano Markidis. Leveraging hpc profiling and tracing tools to understand the performance of particle-in-cell monte carlo simulations. In European Conference on Parallel Processing, pages 123–134. Springer, 2023. [46] Weihang Xu, Maryam Fazel, and Simon S Du. Toward Global Convergence of Gradient EM for Over-Parameterized Gaussian Mixture Models. arXiv preprint arXiv:2407.00490, 2024.