Plasma-PEPSC D1.13 Updated Assessment of Developments and Revision of Scientific Challenge, Predicting Near-Earth Space Dynamics with Vlasiator

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 D1.13 Updated assessment of developments and revision of scientific challenge - Predicting Near-Earth Space Dynamics with Vlasiator WP1: Plasma Simulations - Codes and Grand Challenges Date of preparation (latest version): 30/06/2025 Copyright©2023 – 2027 The Plasma-PEPSC Consortium D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 2 DOCUMENT INFORMATION Deliverable Number D1.13 Deliverable Name Updated assessment of developments and revision of scientific challenge - Predicting Near-Earth Space Dynamics with Vlasiator Due Date 30/06/2025 Deliverable lead UoH Authors Urs Ganse (UoH), Markus Battarbee (UoH), Leo Kotipalo (UoH) and Markku Alho (UoH) Responsible Author Urs Ganse (UoH) E-mail: [email protected] Keywords Plasma simulation, Vlasov equation, mesh refinement, sparse grids, profiling, performance extrapolation WP/Task WP1/Task D1.13 Nature R Dissemination Level PU Final Version Date 30/06/2025 Reviewed by Vassilis Papaefstathiou (FORTH) & Valentin Seitz (BSC) D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 3 DOCUMENT HISTORY Partner Date Comment Version UoH 23/05/2025 First draft, introduction and code description 0.1 UoH 28/05/2025 Development summary and updated gap analysis 0.2 UoH 02/06/2025 Roadmap, summaries and internal proofreading 0.2.1 UoH 03/06/2025 Introduction and conclusion 0.2.2 UoH 03/06/2025 Final draft ready for internal review 0.3 KTH 04/06/2025 Final draft updated for internal review 0.4 UoH 17/06/2025 Revised draft after internal review 0.5 KTH 23/06/2023 Final cleanup for submission 1.0 D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 4 Executive Summary In this deliverable, we provide an update of the progress for the Vlasiator plasma simulation code, in terms of its development progress towards the scientific target case within Plasma-PEPSC. It gives an outline of the basic physical environment and state of the art of near-Earth plasma simulations that Vlasiator performs, and reiterates the scientific target case which would be desirable to simulate within the timeframe of this EuroHPC project. Performance profiling results of Vlasiator have been updated at various scales, and their results are being presented, in direct comparison to the earlier results from deliverable D1.4 [1]. Based on the profiling data, performance extrapolation to the scientific target case shows that Vlasiator currently has a performance gap of a factor 32 compared to what would be required, down from the initial gap value of 52 at the beginning of the project. We discuss further steps in the development of the simulation code to address this gap, and give an overview of the planned improvements to the code. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 5 Contents 1 Introduction 6 2 Background and Scientific Challenge 6 3 Summary of developments 8 3.1 GPUporting ................................. 8 3.1.1 Memory management and architectural improvements . . . . . . . 9 3.1.2 Vlasov solver porting . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.3 Blockadjustment........................... 10 3.1.4 Initialisation and data reductions . . . . . . . . . . . . . . . . . . 10 3.1.5 Unified codebase for CUDA and HIP/ROCm . . . . . . . . . . . . 11 3.1.6 Fieldsolver .............................. 11 3.1.7 Ionosphere solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Non-global timestepping . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3 Load balancing advancements . . . . . . . . . . . . . . . . . . . . . . . . 12 3.4 MPI optimisations based on benchmarking insights . . . . . . . . . . . . 13 4 Results: Revision of Gap analysis 14 4.1 Roadmap ................................... 16 5 Related Work 17 6 Conclusion 17 D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 6 1 Introduction Vlasiator[2, 3] is a kinetic plasma simulation system designed and focused on simulation of the near-Earth space environment. Its scientific goal is to simulate the plasma processes in and around Earth’s Magnetosphere on ion-kinetic scales, to serve both as a fundamental physics laboratory and work in tandem with spacecraft instrument design and data analysis efforts from national and international space agencies. 2 Background and Scientific Challenge Figure 1: Example of a high-resolution magnetosphere simulation output from Vlasiator: Earth, in the centre, generates its magnetic dipole field (field lines are shown in purple), which gets deformed and shaped by the arriving solar wind from the sun (off the right in this figure). The plotted colourmap shows the plasma temperature, as it gets heated by the bow shock and kinetic plasma phenomena. The specific scientific challenge that the Vlasiator team is working towards is a complete, three-dimensional modeling run of Earth’s entire magnetosphere at resolutions that correspond to ion kinetic scales. In Deliverable D1.4 [1] at the start of the PlasmaPEPSC project, we formulated the concrete requirements for our target case, quoted here verbatim from that document: •A simulation of global coupling of Solar Wind↔Magnetosphere↔Ionosphere in a simulation box that contains part of the solar wind and extends to at least −100 RE in the tail. The overall physical size of the simulation box directly correlates with the computational intensity of a simulation. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 7 •The inner simulation boundary should be located at 2.5 RE(or further in). As the magnetic field strength of Earth’s dipole field rises sharply close to Earth (B∼r−3), this creates a numerical challenge: The plasma’s Larmor frequency Ω = qiB mineeds to be resolved by the simulation timestep. Thus, this parameter affects computational scaling by a factor of 1/x3. •Spatial resolution of the simulation mesh (at the highest refinement level) should be 500 km or finer. This means one additional refinement level compared to the current top-end Vlasiator simulations, and a factor of 8 more field solver cells. •The run should continue for at least two simulated hours, about four times as long as the current record. •As opposed to the current homogeneous solar wind inflow conditions, the inflowing solar wind conditions should be time-varying (representing realistic solar wind transient effects, like turbulent CME sheath regions or flux ropes). •The outer size of the simulation domain should be chosen so that it includes Earth’s wave foreshock structure, with its spatially varying and complex structure of velocity distribution functions. •We want to be able to automatically make decisions about dynamic mesh refinement and subgrid modeling at run-time, in order to track the physical phenomena of interest with maximum possible fidelity. At the start of the Plasma-PEPSC project, we identified that all required simulation components and solvers are in place: A suitable ionosphere model [4], adaptive mesh refinement of the spatial grid [3], an appropriately low-diffusivity semilagrangian Vlasov solver [5]. However, the interplay of these components was not fine-tuned, and overall performance made the computational requirements for this scientific goal seem out of reach. The initial gap analysis per D1.4 was determined on the assumption that a realistic large computing grant on an exascale machine (corresponding to half of the machine’s resources for a month) should be sufficient to run the target case. At that point, this corresponded to 800 nodes of LUMI-C. The performance gap for Vlasiator is thus calculated by estimating the extrapolated computational cost of the science case on 800 nodes of LUMI-C, and dividing it by 80 million core-hours: GAP = (∆tVlasov + ∆tFields)×ntimesteps 80 ×106core-hours .(1) In D1.4, the target case’s numerical properties were specified in terms of concrete numbers of simulation spatial grid cells, fieldsolver grid cells, velocity distribution function (VDF) cells and total simulated physical time. In Table 1, we compare these numbers with the current maximum values actually obtained in simulation runs. Note that these are from distinct runs: The longest simulation currently conducted had a much lower spatial resolution, and hence cell count, than the actual high-fidelity science runs. Together, these runs have demonstrated that Vlasiator solvers have the required numerical stability and resolution capability for performing the target case science run. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 8 Quantity Target case Recent prod. runs Fulfilled? Spatial grid cells 5.56 ×1072.26 ×107Almost Fieldsolver cells 4.48 ×1094.71 ×108No VDF cells 4.61 ×1012 4.69 ×1012 Yes Physical run time >2 hours 2h 46min Yes Table 1: Comparison of Vlasiator target case numerical quantities as determined in D1.4 with most recent Vlasiator run results. 3 Summary of developments Vlasiator development in the last two years focused primarily on three large refactoring approaches to gain effective performance on current HPC architectures: Porting of Vlasiator to GPUs (with complete data residency in GPU memory), removal of the global timestep constraints and new algorithms for load balancing. Each of these are outlined below. The fully run-time adaptive mesh refinement of Vlasiator, published in [6], has now been taken into use in production. Maximum resolution is now adaptively following the physics of the simulation, and coarser cells are employed in spots with less scientific interest, freeing up resources dynamically. While this does not effect any measured change in the performance numbers as a function of cell count, like the FOM calculations stipulate, it does provide a significant step towards the target science case, as the required numbers of Vlasov grid cells are significantly reduced for the same physical fidelity. Vlasiator’s code development practices continue to evolve, finding inspiration and support from the other consortium partners. Vlasiator’s primary git repository now performs automatic testing for continuous integration (CI) based on Github actions, with runners hosted on the local cluster of the University of Helsinki (hardware provided courtesy of the Finnish Grid and Cloud Initiative, FGCI) and the HCA cluster at Barcelona Supercomputing Centre. Code review meetings and integration-focused “merge days” have become central to the development progress. In addition, a multitude of small performance micro-optimisations and removal of individual performance hot spots have been performed. 3.1 GPU porting A GPU-compatible port of Vlasiator has been envisioned since the late 2010s, with preliminary work towards directive-based OpenACC porting taking place in 2019 and 2020. A dedicated push towards porting critical functionality of Vlasiator to CUDA and HIP/ROCm has been underway since the beginning of 2021. Due to the object-oriented programming approach and non-uniform non-constant grid structures of Vlasiator and the way the physical solvers interact, much effort was required in preparatory and groundbreaking work, in order to facilitate the required operations to take place on GPU hardware. In particular during the first 30 months of Plasma-PEPSC, this work has come to fruition. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 9 3.1.1 Memory management and architectural improvements Vlasiator utilizes a sparse velocity space [5] where for each spatial cell, the size and shape of the modeled and stored region of phase-space is adjusted several times per time step. This is crucial for being able to fit the physical problem into memory, with sparsity accounting for over 99% of savings in memory. As such, when acting upon a spatial cell (represented by an object with various data containers), the sizes of these containers can vary from cell to cell and from timestep to timestep. Consequently, information of these cell objects needs to be retained both on-device (for the computational kernels to be able to adjust and update the relevant information) and on the host (for the program to be able to define suitable launch parameters for computational kernels). A crucial piece to the puzzle of being able to manage these containers both on host and on device was the in-house development of the portable unified memory vector and hash map constructs [7]. These are built to support both CUDA and HIP/ROCm and allow the correct data to be accessible at all times from both host and device sides. During the winter of 2024–2025, an approach was finalised which allows crucial data to remain in unified memory but reside as only paged to device during regular operations. Host-side copies or cached counters are used for kernel launch parameters, and these values are carefully maintained as up-to-date. As the data still resides in unified memory, it is trivial to perform sanity checks between cached and actual values. This caching approach was crucial in preventing the memory subsystem from paging this unified memory data between the host and the device. Paging of data back and forth has a significant detrimental effect on performance. 3.1.2 Vlasov solver porting As of spring 2025, the computationally intensive portions of the Vlasov solvers of Vlasiator act on data residing in either device memory or unified memory which does not need paging to host. As such, the basic performance criteria for on-device computation are fulfilled. The translation solver (spatial propagation) has been optimised and kernel merging has been performed, so that for each Cartesian translation direction there is only a single kernel launch which processes the whole necessary data, utilising pre-allocated temporary buffers for the Vlasov remapping. There are still some remaining optimisations pending, such as implementing more fused multiply-add (FMA) instructions. Also, the current memory access patterns suffer from significant long scoreboard stalls, a topic for future investigation and optimisation, as the reason for this is not yet known. The translation solver is, however, well-optimised in comparison with other portions, and has been in a relatively stable state. The semi-Lagrangian acceleration solver (based on the SLICE-3D approach of [8]) was ported to a CUDA kernel as early as January 2023. Many of the auxiliary tasks of preparing data for acceleration were ported piece-by-piece, but the legacy of serial CPU operations was still present for a long time. During the spring of 2025, the rest of the acceleration solver has been re-designed with more compact kernels which utilise best practices of GPU programming and will facilitate increased parallelism. This re-writing is still underway at the time of this report. A crucial step in this task was actually re-designing the algorithms for parallel operations. For example, instead of utilising a CPU or GPU D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 16 Figure 6: Scaling graph and extrapolation of Vlasiator’s figure of merit (FOM) calculated according to the calculation defined in the grant agreement (eq. 2). The data points and colours are identical to Figure 5. Additionally, the success criterion outlined in the GA (FOM <1000) is drawn as a dashed horizontal line. 4.1 Roadmap The next major steps in Vlasiator development consist of taking into production the recent advances outlined in section 3 above. GPU support has been merged into the primary development branch, but appropriate resources for a GPU-based production run have not yet been obtained, as the scaling behaviour was still varying strongly during development. This will now be remedied. With proper GPU scaling data, quantitative decisions can be made whether to target the next production runs to GPU or CPU systems. The non-global timestep implementation, is still in a prototype stage. It will require thorough testing and integration runs, before its performance impact at scale can be quantified. In any case, it will be a substantial step towards moving the simulation’s inner boundary further earthward. This process is expected to be completed throughout the year 2026. Yet, the remaining gap value of 31.87 will not be sufficiently reduced by these two developments alone. As the scaling analysis in Section 4 demonstrates, the field solver presents a scaling barrier. Work on GPU porting of this component of the code is well underway (compare Section 3.1.6) and an advanced MPI partitioning scheme, in which only a subset of the participating compute nodes solve the Fields is already functional, but has not been tested at scale. In this problem, we are in exchange of ideas with the GENE and GENE-X development teams, which likewise focus on field solver scalability progress. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 17 5 Related Work As part of the collaborative effort that is the Plasma-PEPSC project, Vlasiator development happens in direct correspondence and interaction with the other consortium partners through the technical work packages (WP2–WP5). Consequently, the deliverables associated with these workpackages’ tasks contain more detailed technical information about the individual performance optimisation, benchmarking and consolidation work that has been performed. Scientific analysis of the production-scale magnetospheric runs that were used in the scaling analysis is still ongoing. First publications based on their data have recently been published in [20, 21, 22]. 6 Conclusion Vlasiator’s performance development towards the exascale target case progresses as the Plasma-PEPSC project continues. The performance gap for the scientific target case has shrunk from 52 at the beginning of the progress to a new value of GAP = 31.87 as of June 2025. Several large code refactoring steps have been completed, most prevalently efficient GPU support, non-global timestepping and significantly improved load balancing. These are being integrated, tested and benchmarked for production readiness. Vlasiator has shown numerical solver stability for the required physical time duration and resolution of the target case. Further developments are already ongoing, and the development appears to be on track to fulfill the target science case by the conclusion of the project. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 18 References [1] Urs Ganse. Plasma-pepsc d1.4 initial definition of the scientific challenges and initial evaluation plan - predicting near-earth space dynamics with vlasiator, June 2023. [2] Minna Palmroth, Urs Ganse, Yann Pfau-Kempf, Markus Battarbee, Lucile Turc, Thiago Brito, Maxime Grandin, Sanni Hoilijoki, Arto Sandroos, and Sebastian von Alfthan. Vlasov methods in space physics and astrophysics. Living Reviews in Computational Astrophysics, 4(1):1, August 2018. [3] Urs Ganse, Tuomas Koskela, Markus Battarbee, Yann Pfau-Kempf, Konstantinos Papadakis, Markku Alho, Maarja Bussov, Giulia Cozzani, Maxime Dubart, Harriet George, Evgeny Gordeev, Maxime Grandin, Konstantinos Horaites, Jonas Suni, Vertti Tarvus, Fasil Tesema Kebede, Lucile Turc, Hongyang Zhou, and Minna Palmroth. Enabling technology for global 3d + 3v hybrid-vlasov simulations of near-earth space. Physics of Plasmas, 30(4):042902, April 2023. [4] Urs Ganse, Yann Pfau-Kempf, Hongyang Zhou, Liisa Juusola, Abiyot. Workayehu, Fasil Kebede, Konstantinos Papadakis, Maxime Grandin, Markku Alho, Markus Battarbee, Maxime Dubart, Leo Kotipalo, Arnaud Lalag¨ue, Jonas Suni, Konstantinos Horaites, and Minna Palmroth. The vlasiator 5.2 ionosphere – coupling a magnetospheric hybrid-vlasov simulation with a height-integrated ionosphere model. Geoscientific Model Development, 18(2):511–527, 2025. [5] S. von Alfthan, D. Pokhotelov, Y. Kempf, S. Hoilijoki, I. Honkonen, A. Sandroos, and M. Palmroth. Vlasiator: First global hybrid-vlasov simulations of earth’s foreshock and magnetosheath. Journal of Atmospheric and Solar-Terrestrial Physics, 120(0):24 – 35, 2014. [6] Leo Kotipalo, Markus Battarbee, Yann Pfau-Kempf, and Minna Palmroth. Physicsmotivated cell-octree adaptive mesh refinement in the vlasiator 5.3 global hybridvlasov code. Geoscientific Model Development, 17(16):6401–6413, 2024. [7] K. Papadakis, M. Battarbee, U. Ganse, Y. Pfau-Kempf, and M. Palmroth. Hashinator: A portable hybrid hashmap designed for heterogeneous high performance computing. Frontiers in Computer Science, 6, 2024. [8] M. Zerroukat and T. Allen. A three-dimensional monotone and conservative semilagrangian scheme (SLICE-3D) for transport problems. Quarterly Journal of the Royal Meteorological Society, 138(667):1640–1651, 2012. [9] 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. Journal of Physics: Conference Series, 2997(1):012010, April 2025. [10] Ga¨el Guennebaud, Benoˆıt Jacob, et al. Eigen v3. http://eigen.tuxfamily.org, 2010. D1.13: Predicting Near-Earth Space Dynamics with Vlasiator 19 [11] E. G. Boman, U. V. Catalyurek, C. Chevalier, and K. D. Devine. The Zoltan and Isorropia parallel toolkits for combinatorial scientific computing: Partitioning, ordering, and coloring. Scientific Programming, 20(2):129–150, 2012. [12] Martin Schulz, Jeremy Williams, Markus Battarbee, Leo Kotipalo, Felix Jung, Srikanth Sathyanarayana, Isaias Urena, and Tilman Dannert. Plasma-pepsc d3.2 updated report on mpi, load balancing and resilience in plasma codes, December 2024. [13] Herman Haverkort. How many three-dimensional hilbert curves are there? Journal of Computational Geometry, 8(1):206–281, 9 2017. [14] Leo Kotipalo, Markus Battarbee, Yann Pfau-Kempf, Vertti Tarvus, and Minna Palmroth. Load balancing in strongly inhomogeneous simulations – a vlasiator case study, 2025. https://arxiv.org/abs/2505.20908. [15] Michael Wagner, Stephan Mohr, Judit Gim´enez, and Jes´us Labarta. A structured approach to performance analysis. In Christoph Niethammer, Michael M. Resch, Wolfgang E. Nagel, Holger Brunst, and Hartmut Mix, editors, Tools for High Performance Computing 2017, pages 1–15, Cham, 2019. Springer International Publishing. [16] Extrae Team. Extrae |BSC-Tools. [17] V. Pillet, J. Labarta, A. Cortes, and S. Girona. Paraver: A tool to visualize and analyze parallel code. In World occam and Transputer User Group Technical Meeting, page 17, 04 1995. [18] Valentin Seitz, Markus Battarbee, Urs Ganse, Vertti Tarvus, Minna Palmroth, and Marta Garcia-Gasulla. Performance analysis of vlasiator at production-scale. In Proceedings of the Third EuroHPC user day, 2025. submitted. [19] Sebastian von Alfthan. Phiprof – parallel hierarchical profiler. https://github. com/fmihpc/phiprof, 2023. [20] Lauri P¨ank¨al¨ainen, Giulia Cozzani, Markus Battarbee, Urs Ganse, Yann PfauKempf, Jonas Suni, and Minna Palmroth. Identifying magnetotail jet front signatures in a 3d + 3v global hybrid-vlasov simulation. Journal of Geophysical Research: Space Physics, 130(6), May 2025. [21] Giulia Cozzani, Markku Alho, Ivan Zaitsev, Hongyang Zhou, Sanni Hoilijoki, Lucile Turc, Maxime Grandin, Kosta Horaites, Markus Battarbee, Yann Pfau-Kempf, Urs Ganse, Kostis Papadakis, and Minna Palmroth. Interplay of magnetic reconnection and current sheet kink instability in the earth’s magnetotail. Geophysical Research Letters, 52(2), January 2025. [22] I. Zaitsev, G. Cozzani, M. Alho, K. Horaites, H. Zhou, A. Kit, Y. Pfau-Kempf, S. Hoilijoki, U. Ganse, M. Battarbee, K. Papadakis, J. Suni, M. Dubart, F. TesemaKebede, A. Workayehu, V. Tarvus, L. Kotipalo, V. Koikkalainen, L. Turc, and M. Palmroth. Ion-mediated tearing and kink instabilities in the earth’s magnetosphere: Hybrid-vlasov simulations. Journal of Geophysical Research: Space Physics, 130(1), January 2025.