Porting the grid-based 3D+3V hybrid-Vlasov kinetic plasma simulation Vlasiator to heterogeneous GPU architectures (2024)

1.1 Vlasiator

Vlasiator utilizes explicit solvers with a Strang splitting [Strang, 1968, Einkemmer and Ostermann, 2014] approach, decomposing evolution of (1) into six Cartesian shear operations; three in velocity space (referred to as acceleration), and three in spatial space (referred to as translation). The Lorentz force effect on the plasma distribution (first and third terms of equation 1) can be (in the non-relativistic case) described as gyration in a frame of reference where the motional electric field cancels out the existing electric field. Using the SLICE-3D [Zerroukat and Allen, 2012] algorithm, this is deconstructed into three Cartesian direction aligned shear transformations.

1.2 Discretization

Vlasiator discretizes ion VDFs on Cartesian grids, storing and propagating phase-space average values in units of $(\Delta x\Delta y\Delta z\Delta v_{x}\Delta v_{y}\Delta v_{z})^{-1}$. The simulation domain consists of a 2D or 3D spatial Cartesian mesh, which can be either uniform, or use cell-based mesh refinement with an octree approach. For 2D simulations, the third dimension consists usually of a thickness of a single cell with periodic boundary conditions. Non-unitary thicknesses (so-called 2.9D simulations with a cylindrical inner boundary) are also supported, such as those used for a study of dayside reconnection with multiple X-lines [Pfau-Kempf etal., 2020]. When simulating in three spatial dimensions, the octree approach is used to refine low-resolution cubic cells into 8 children each with half the dimensions of the parent. Within a neighbourhood of 3 cells in each direction, refinement level can change by only one level, to prevent too drastic resolution boundaries. Refinement can be repeated several times, with current simulations utilizing up to 4 different resolution levels. Details of the AMR approach used in Vlasiator are presented in Papadakis etal. [2022], Ganse etal. [2023] and Kotipalo etal. [2024]. Typical Vlasiator 2D (or 5D; 2D+3V=5D) simulation spatial resolutions have been between $227$ and $300\,\mathrm{km}$, slightly above the ion inertial length in the solar wind, but capable of modelling ion-kinetic effects as shown in Pfau-Kempf etal. [2018]. In 6D runs, the highest spatial resolution in regions of interest has been $1000\,\mathrm{km}$ which is above the ion inertial length, but those simulations have been shown to still include significant ion-kinetic dynamics [Grandin etal., 2023, Palmroth etal., 2023].

Each spatial cell on the spatial mesh contains a local representation of phase-space, also discretized on a Cartesian velocity mesh. To conserve memory and computational resources, Vlasiator utilizes a sparse velocity space approach, only storing and propagating those regions of velocity space which contain non-negligible phase-space densities. This is performed by declaring the maximum allowable velocity space extent and dividing it into cubic blocks of $(4\times 4\times 4)$ or, as implemented recently, $(8\times 8\times 8)$ velocity space cells each. Internally this is declared with notation of WID=4 or WID=8, with values of WID3=64 and WID3=512, respectively. Typical maximal velocity space extents have been $\pm 4.02\cdot 10^{6}\,\mathrm{ms}^{-1}$ in each direction with 67 blocks per direction, resulting in a velocity space resolution of $30\,\mathrm{kms}^{-1}$, or with 6D simulations, $\pm 4\cdot 10^{6}\,\mathrm{ms}^{-1}$ with 50 blocks per direction, leading to a velocity space resolution of $40\,\mathrm{kms}^{-1}$. This reduction of velocity resolution hand-in-hand with spatial resolution follows the findings of Pfau-Kempf etal. [2018]. Active blocks are added to the velocity space or removed from it as necessary, as described in von Alfthan etal. [2014]. We ephasize that the term sparse is used in this manuscript to discuss this mechanism for partitioning and selectively storing and propagating portions of potential velocity space, and thus it is not directly related to e.g.sparse matrix algorithms. An example of a Vlasiator simulation domain, showcasing spatial mesh refinement and differing spatial cell computational weights due to varying block counts, is shown in Figure 1.

2.1 Sparse block representation in memory

The VDF at a given spatial location is represented by two C++ classes: a Velocity mesh (Vmesh) and a Velocity Block Container (VBC). The Vmesh links memory buffer positions (local address offset i.e.the local block identifier, LID) with physical velocity coordinates (global block identifier, GID). It contains a hash map (globalToLocalMap) with GIDs as keys and LIDs as values. The complementary data structure is a vector (localToGlobalMap) with size equal to that of the map, and where each index (LID) stores the corresponding GID. The VBC, on the other hand, contains the actual VDF phase-space data $f_{s}$ (Block data container, vector) as WID3 floating point phase-space density values for each block and LID, as well as pre-calculated metadata (Block parameters, vector) for each WID3-sized block in question. LIDs are thus simply indices into WID3-staggered memory space, whereas GIDs represent a position in velocity space, with values starting from 0 for the smallest velocity value in all Cartesian directions, incrementing first along the $V_{x}$ axis, then along the $V_{y}$ axis, and finally along the $V_{z}$ axis. When evaluating 6-dimensional phase-space neighbours of any given block, spatial neighbours will have the same GID, whereas velocity neighbours will have a GID differing by $\pm 1$, $\pm N_{\mathrm{b},x}$, or $\pm N_{\mathrm{b},x}N_{\mathrm{b},y}$ where $N_{\mathrm{b},i}$ is the maximum possible number of blocks in Cartesian direction $\mathrm{i}$.

2.2 Block adjustment requirements

A naïve estimation of memory consumption for a 6D Vlasiator simulation with fully populated phase-space grids is $3.5\times 10^{15}$ phase space cells, a memory requirement of about 14 petabytes [Ganse etal., 2023, Kotipalo etal., 2024]. A significant contributor to Vlasiator’s success is the sparse velocity mesh, reducing the computational domain to include only those regions of phase-space where non-negligible density $f_{s}>f_{\mathrm{th}}$ is found, yielding a memory reduction of over two orders of magnitude.This strength comes with the cost of having to constantly adjust and update the structure of phase-space stored in each cell. The acceleration solver evaluates future velocity space extents after performing each shear operation, flagging content blocks which are no longer needed and those which need to be created. The translation solver on the other hand requires that neighbouring spatial cells contain the required data structures for phase-space density to flow into. These requirements are described as the ghost halo regions, as indicated in Figure 3. After running the acceleration solver, each content block is flagged as being either above the sparsity threshold $f_{\mathrm{th}}$ (i.e.having significant content), or below it (i.e.not having significant content). Then, the halo requirements are evaluated for all content blocks, both for velocity and spatial neighbours. Finally, required halo blocks are created, and non-content blocks not required by any halos are deleted.

Managing velocity space sparsity via blocks allows Vlasiator to use standard data units of size WID3, useful for vectorization and GPU approaches. Vlasiator has now been amended to support blocks of size $(8\times 8\times 8)$ (WID3=512) in addition to the longstanding $(4\times 4\times 4)$ (WID3=64) block size. This results in much less overhead from creating and deleting blocks since a smaller number of blocks is required to describe any given phase-space distribution function, assuming the resolution remains the same. At the same time, the halo requirements will, on average, lead to a larger memory footprint and more actual phase-space cells to be propagated. This is a trade-off, but one which is expected to prove beneficial on GPU architectures. The implications on halos of changing the value of WID for a simple configuration, whilst maintaining constant maximum velocity space extents, is exemplified in Figure 4.

2.3 Other grid objects in memory

In addition to the DCCRG grid containing spatial cells with particle species VDF data, each MPI task also owns several FSgrid objects which house uniform Cartesian rectangular domains of scalars or vectors for field solver use. The total spatial domain is decomposed uniformly over all MPI tasks, with addition ghost buffer regions around local data, which can be updated from the up to 26 neighbouring tasks. The division into several different FSgrid objects is to facilitate ghost communication of only the necessary data for each step of the field solver operation. The Vlasiator field solver, consisting of much smaller data domains than the Vlasov solver, has yet to be ported to heterogeneous GPU architectures.

A final data object of Vlasiator is the ionospheric grid [Ganse etal., 2024], which is contained on every task in its entirety and does not contain any MPI interface. The ionospheric grid exists only when simulating a magnetospheric run with an ionospheric inner boundary, and its solvers have yet to be ported to GPU architectures.

2.4 Heterogeneous memory decomposition of data structures

As Figure 2 introduced in Section 2 shows, the data structures of Vlasiator are non-trivial when considering heterogeneous memory architectures such as those encountered in modern GPU-powered supercomputing nodes [Zheng etal., 2016, Wang etal., 2014]. Computational activities performed by the CPU can usually only access host-side memory, and parallel computation kernels utilizing GPU hardware can only efficiently access device-side memory. Current pre-exascale and exascale supercomputers are built upon node architectures which house discrete GPUs controlled by standard CPUs, each with discrete memory domains.GPUs are built using high-bandwidth memory [HBM, Zhu etal., 2018], allowing rapid availability of data for GPU computations with speeds of up to 3.2 TB/s on MI250x and 2 TB/S on NVIDIA A100. A bottleneck, however, is the connection between GPU and CPU hardware, currently limited tospeeds of tens to a few hundreds of GB/s per direction (depending on NVLink or PCI Express implementation),necessitating planned and controlled movement of computation-critical data between the two memory subsystems. Upcoming hardware such as the Nvidia Grace Hopper GH200 will upgrade this to a 900 GB/s bidirectional NVLink-C2C connection as well as a transparent unified memory domain doing away with the need for data migration. GPU HBM bandwidth will still remain superior at 3.0 TB/s. Accelerated processing unit (APU) chips such as the AMD Instinct MI300A will contain a single memory domain [Vijayaraghavan etal., 2017], which may simplify memory management on future exascale architectures whilst simultaneously providing massive bandwidth (5.3 TB/s on the MI300A).

GPU application programming interfaces (APIs) provide access methods allowing various ways of copying data between host and device memory regions. A user-friendly implementation (supported on modern architectures such as NVIDIA Pascal and later and AMD MI200+, which still have physically distinct memory domains) is that of Managed Memory or Unified Memory (UM), where a single virtual memory address can point, as necessary, to either host or device memory. The UM subsystem will utilize a paging mechanism, and can even allow oversubscription of device memory, paging (migrating) memory to and from the device as necessary. Whenever the CPU or GPU needs to access data residing in UM, the paging system will be queried for the locality of the data. If the data is not in the required region, a page fault occurs, triggering migration of the data to the correct domain. To reduce page faulting (which can harm performance), users can prefetch data manually to the desired domains. One noteworthy aspect of prefetching and page faulting is that the UM system only migrates data in units of the page size, which can lead to unintended migration of data units residing within the same page as an intentionally migrated data unit. However, the benefit of unified memory is that an unit of data which is accessed from both host and device is ensured to always contain the correct value, irrespective of whether the user has declared prefetches. Should a data construct contain separate host-only and device-only memory buffers, it is the responsibility of the user to perform manual memory copies in order to maintain correct values for host and device activities, respectively. Additionally, automatic oversubscription of device memory is not possible if separate buffers are used. When porting a large simulation software such as Vlasiator, this reliability of UM was deemed a high priority, and thus, many Vlasiator data structures were migrated to UM instead of maintaining pure host and device copies of data buffers.

Of the data containers described above and showcased in Figure 2, only the lowest level of data containers, i.e.the hashmaps and vectors (dark and light blue in the figure) reside in unified memory. This includes both the data buffers and also metadata such as the container sizes. For FSgrid data, both unified memory and dual host-and-device memory versions have been implemented. To facilitate access to Vlasov data structures within the Vmesh and VBC objects from both host and device environments, Vlasiator stores discrete host and device copies of the actual Vmesh and VBC objects, which both point to the same hash maps and vectors residing in UM.

2.5 Data structure re-ordering for Vlasov solver use

The Vlasiator Vlasov solver approach, introduced in Section 1.1,is to decompose the Vlasov equation into six one-dimensional Cartesian shear or advection operations. Each advection operation takes a number of sampling points in the relevant direction, using them to fit a polynomial which describes the variation of phase-space density $f_{s}$ along the direction of propagation. This polynomial is then advected with relevant speed and direction and integrated over the target phase-space cell dimensions.

For each advection thus performed by our semi-Lagrangian approach, we require a varying number of control points depending on the polynomial reconstruction order. In order to limit the number of ghost cells required, we usually constrain translation to use parabolic polynomials for each cell, whereas acceleration utilizes quadratic polynomials. This results in each advection calculation requiring both the local phase-space density $f_{s}$ value and values of a number of neighbours in the direction of advection. For this reason, the first step of each Vlasov solver cycle is to take the existing Vlasov data structure andre-arrange it so that memory access operations will occur on contiguous and squential data regions.

For spatial translation, the solver will act on one global block identifier (GID) at a time, collecting the data along the translated direction, utilizing the pencil formalism described in Ganse etal. [2023]. For each spatial cell, the GID being translated will be linked to the local block identifier (LID) with the globalToLocalMap hashmap object inside the Vmesh. The block data from each memory address given by the LID is then copied into a temporary device-only buffer, providing the solver with aligned block data along the direction of propagation.The original data buffers pointed to by each LID are then emptied (written to zero), and the semi-Lagrangian propagator is then able to write its output into the original data locations.

For velocity-space acceleration, phase-space densities will be advected along one Cartesian velocity direction at a time. The GID associated with each data block is such that values increase in the order $V_{x}\rightarrow V_{y}\rightarrow V_{z}$. In order to ensure the block data is aligned along the direction of advection $i$, we compute a transposed $\widetilde{GID}_{i}$ where $\widetilde{GID}_{x}=GID$. $\widetilde{GID}_{y}$ is defined by the order $V_{y}\rightarrow V_{z}\rightarrow V_{x}$ and $\widetilde{GID}_{z}$ by the order $V_{z}\rightarrow V_{x}\rightarrow V_{y}$. In order to build the ordered column data structures, we utilize a sorting algorithm on the $\widetilde{GID}$-values, thus gaining access to the correctly sorted LIDs as well. The data from these addresses is then copied to a temporary buffer, which can be fed to the Vlasov solver. In the meanwhile, the contents of the velocity blocks are reset, so that the Vlasov solver can write the results of the solver operation into the cell data contents.

3.1 Directive-based approaches

Vlasiator, utilizing existing CPU-based computational facilities to a great degree, is built to perform calculations both in parallel distributed over MPI to several tasks, as well as running several OpenMP threads in parallel, possibly also benefiting from hyperthreading. As such, OpenMP directives are implemented in several parts of the Vlasiator codebase, and directives were an interesting direction for trialling GPU porting. During the autumns of 2019 and 2020, the Vlasiator development team participated in Hackathons where OpenACC directive-based GPU programming was tutored and implemented. However, our testing suggested that OpenACC was not going to provide us with sufficient throughput, as mere pragma-based offloading of selected code sections would not easily result in sufficient re-factoring or provide us with the granularity of memory management foreseen to be required. Vlasiator, implementing a 6-dimensional advection equation, requires each point of VDF data to be available for each solver step, save for certain ghost cell information, meaning VDF data must remain on the device. Additionally, compiler support of pragma-based approaches during our initial porting efforts was lacking. As such, we did not pursue OpenACC approaches further, nor did we trial OpenMP 5.0 GPU-based directives.

3.2 Abstraction layers

A useful approach in porting existing software to new architectures is to utilize a portable abstraction layer, which can be externally maintained and updated to allow interfacing with new hardware and programming languages without the end-user needing to change their code. Abstraction can be built into the grid library itself [such as in AMReX, Zhang etal., 2019] or provided through an external portability layer, which in turn may depend on libraries for memory management. Examples of such portability layers include RAJA (Beckingsale etal. 2019 with array abstractions via CHAI and memory management via UMPIRE Beckingsale etal. [2020]), alpaka (Matthes etal. 2017 with memory management via LLAMA) and KOKKOS (Trott etal. 2022 consisting of both core and tool packages). Industry giants are also attempting to tackle the issues of vendor lock-in: for example the Intel oneAPI [Costanzo etal., 2021, Christgau and Steinke, 2020] implements specialized APIs as well as the portable SYCL language [Alpay etal., 2022]. An interesting note is also that HIP which is often viewed as the de facto language on AMD hardware, is in fact an abstraction layer which translates commands to the ROCm interface used by AMD hardware, whilst also supporting translation to CUDA on Nvidia hardware. Larger adoptance of HIP remains to be seen, though efforts of supporting HIP on top of oneAPI do exist [Videau, 2022].

Whilst utilization of abstraction layers brings with it many benefits, it at the same time detaches the developer from the underlying programming language. Algorithms are restricted to constructions supported by the abstraction layer, which may or may not be useful for the task at hand. Due to the inherent complexity of Vlasiator’s memory layout, consisting of a sparse ever-evolving velocity mesh in each spatial cell, and the need to maintain all data in device memory for every solver step, pre-existing abstraction layers were not included in the roadmap of Vlasiator GPU development. After finding directive-based approaches to be too limited, it would have been counterproductive to run into the same issue with portability layers. Abstraction layers are a possible future avenue for fixed-size grids such as the ionospheric and field solver grids, though.

4.1 Unified memory

The Vlasiator simulation code includes several different solvers and routines which act on the modelled particle VDFs. Initialization routines set up the simulation and boundary cell management can re-write cell contents. VDFs are reduced to obtain plasma velocity moments for use by the field solver. Block values are evaluated to find maximum allowed timesteps due to spatial translation. Output routines can process the VDFs to extract partial moments or precipitating energy spectra. Due to all of these accessors, the VDF is stored in unified memory, which can be paged to host or device as necessary. Using unified memory may incur a performance penalty, but reduces the risk of one code section acting on data which might be out-of-date (having been altered by another part of the code) or trying to access inaccessible memory addresses. During early development, small subsections of the solvers were ported one-by-one, allowing verification of correct code behaviour without having to convert all functions to GPU kernels at once. As data movement between the host and the device is much slower than the HBM of GPU hardware, and would cause a significant performance bottleneck, the final state of the Vlasiator code has all VDF-accessing routines run through GPU kernels.

One feature of unified memory on current NVIDIA and AMD architectures is the possibility of page faulting data to the required location. Unified memory is managed as pages, large units of data (e.g.4 kilobytes), which reside either in host or device memory. When accessing such managed memory, a lookup to the UM page table is performed, and assuming the data is already in the correct place, it is read from the address indicated by the page. If the memory resides on the host when requested by the device, or vice versa, a page fault occurs, and the page lookup table access triggers a memory copy, migrating the page to the requested memory domain. The lookup table is then updated to reflect this new status. Page faults are what allow unified memory to be accessed from both host and device, but repeated page faults cause significant slowdowns as computation must wait for the page transfer to complete.

Memory can be prefetched by the user if one knows that the data resides in the incorrect location, however, this incurs an overhead of reading the page lookup table and returning only once the potential transfer of each page associated with the buffer in question has been requested. To completely avoid page faults, a synchronization call requesting for this transfer to complete is also required. The combination of these lookup table accesses and GPU API calls will also incur an overhead, so prefetching should only be performed when the software is sure the data resides in the incorrect memory regime. The current GPU version of Vlasiator contains a number of small counters which reside in unified memory, and the page faulting of these back-and-forth is quite detrimental to performance, but manual prefetching is even worse due to the overhead of GPU API calls. If manual prefetching is performed, the GPU driver may be unable to optimize whether a unit of data should be paged back and forth, or if it should simply be copied into a temporary buffer, maintaining the page location where it is.

Of course if a unit of data is not actually required on both host and device, use of unified memory should be avoided. Pure device memory access does not pass through page lookups, whereas host memory usually does. To further improve host performance, in particular when copying data between host and device buffers, the use of pinned memory is crucial. Pinned memory is declared outside the paging system, bypassing the need to go via host page lookups when communicating with the GPU device. In Vlasiator, we use pure device buffers and pinned host memory for any single-location temporary buffers which allow it.

4.2 Re-using device and unified memory buffers

In addition to managing paging of the unified memory subsystem, care must be taken in managing allocations and deallocations. Within the current CUDA and HIP/ROCm ecosystems, managed or unified memory can only be allocated or deallocated when no GPU kernels or memory transfer operations are active. Thus, these are blocking operations which can have a significant detrimental effect on performance. Thus, the use of buffers should be designed as to minimise allocations and deallocations, re-using buffers where necessary. This is of particular importance to Vlasiator, where required buffer sizes change constantly as the velocity meshes grow and shrink. Vlasiator thus utilizes a method of maintaining slightly larger buffers than required, re-using them for subsequent kernels, and whenever a re-allocation takes place, making the new allocation larger by a factor of e.g. 30% in order to postpone the need of future re-allocations. Allocations of pure device memory (or host memory) do not block the whole device, as they can be performed asynchronously in one GPU stream, and thus they are not as crucial to avoid.

4.3 CPU threads and GPU streams

Vlasiator on CPUs is massively parallel, utilizing MPI for communication between different tasks and supercomputer nodes, OpenMP threading for parallel processing on several CPU cores, and vectorized instructions to process several units of data with a single instruction (SIMD). GPU hardware is also capable of running several operations in parallel, both as GPU warps or wavefronts (SIMT, single instruction multiple threads), but also as several different kernels running in parallel. Vlasiator utilizes this capability in a heterogeneous manner by using several CPU threads to prepare data and queue operations, with the GPU performing the heavy lifting. Each CPU OpenMP thread is given an unique GPU stream to which it can queue operations with each stream acting independent of other streams. Should there be a data dependency between different streams, this can be managed by synchronizing the streams in question, or usually, the whole device, before continuing with enqueueing new GPU kernels or operations.

As each GPU stream will contain e.g.the different sub-kernels required by the acceleration Vlasov solver, and any operations queued to that stream will be performed in order, it makes sense to ensure that each stream has its own re-writeable buffers to act on. Thus, on the host side, each CPU thread (controlling a designated GPU stream) maintains its own temporary buffers. These buffers are re-used when that thread (and, by extension, stream) moves on to the next unit of data to be processed, such as the next spatial cell. Due to the use of these per-thread buffers, the translation solver is also launched in several streams, each using the buffer associated with that thread/stream.

4.4 A Sparse velocity mesh on GPU devices

In order to be able to access and edit the sparse velocity meshes used in Vlasiator from on-device kernels, a hashmap supporting on-device operations was required. During the porting of Vlasiator to heterogeneous GPU architectures, we performed porting bit-by-bit, requiring thus a hashmap which could be accessed from both host and device, in order to provide correct functionality even before all operations had been completed. For this purpose, we developed Hashinator, a portable hashmap library using unified memory and supporting both CUDA and HIP/ROCm. In addition, we made available a portable unified memory vector class called Splitvector, on which Hashinator is built. These tools, along with their portability features and excellent performance, are introduced in Papadakis etal. [2024a] and available at Papadakis etal. [2024b].

Utilizing the Hashinator and Splitvector libraries, the Vmesh and VBC contents (hashmaps and vectors in dark and pale blue in Figure 2) were converted to unified memory constructs, capable of being accessed from both host and device. Though these intermediate objects themselves could also reside in unified memory, we found that it was preferred to maintain host and device copies of Vmeshes and VBCs, which contain host and device copies of splitvectors and hashmaps, which in turn contain unified memory arrays and counters. This helps to reduce page faulting between host and device, but requires the device copy of the objects to be updated whenever a re-allocation occurs in order to update unified memory pointers. A common occurrence when launching kernels is to evaluate the size of the velocity mesh, which would usually mean a call to the size() function of the globalToLocalMap hashmap or the localToGlobalMap splitvector. We found that carefully caching this value as a host-side integer and updating it whenever the velocity mesh size is adjusted resulted in reduced page faulting.

4.5 GPU thread parallelism

As Vlasiator already utilizes memory stored in units of WID3 floating-point values, acting on 32 (CUDA warp, NVIDIA) or 64 (ROCm wavefront, AMD) elements in one pass might have been relatively straightforward. However, the acceleration solver in particular acts only on one WID$\times$WID slice at a time, in order to ensure thread-safe storage of phase-space density values. This small slice-by-slice approach is sub-optimal on GPU hardware if the standard value of WID=4 is retained. A future task will be to refactor the memory access patterns to allow a whole WID3 block to be processed at once, but for now, a transition to WID=8 means that each WID$\times$WID slice has 64 entries and can be efficiently processed on GPU hardware. When a GPU kernel performs operations on a globalToLocalMap hashmap, we utilize the warp accessors of Hashinator, doing lookups of a whole warp length at once, further improving thread-parallelism.

4.6 Porting the Vlasov solvers

The actual process of porting Vlasiator’s Vlasov solvers to GPU kernels is relatively straightforward. Care must be taken to minimize memory transfers between host and device. Re-use of existing buffers is preferred over allocations and de-allocations, which can improve or hinder memory management challenges.

The acceleration solver has to evaluate the contents of the local velocity mesh and re-assemble it into contiguous memory domains dependant on the Cartesian direction in which the shear in question is taking place. One subsection of this involves converting the velocity mesh GIDs to direction-dependent values $\widetilde{GID}$(see Section 2.5). The new $\widetilde{GID}$-values are then sorted in order to achieve contiguous input data to the remapping algorithm.Notably, this is a place where existing sorting libraries such as CUB::sort on CUDA and hipCUB::sort on HIP can be utilized. When evaluating extents of contiguous columns of blocks to be passed to the acceleration solver, thread-parallel GPU search operations provide performance improvements.

The translation solver has a relatively simple flow, allowing merging the preparatory, advection, and storage sections into one large kernel which can be called in parallel utilizing several temporary buffers. The acceleration solver, on the other hand, involves much more transposition and restructuring of data, resulting in a number of small kernels called one after another, which impacts performance by accruing API overhead and sometimes requires passing counter values between the host and the device. After preparatory work, the accelerated cell velocity space must be adjusted to accommodate the new positions in velocity space, which requires calls to the block adjustment algorithms. Furthermore, each spatial cell is accelerated independently of all others, meaning that though they can be processed in parallel, each one requires its own kernel calls in order. Each CPU thread is able to process one cell at a time in its own GPU stream, but extracting acceptable performance from modern GPU hardware may require even further re-factoring of the acceleration solver ordering.

5.1 Serial CPU block adjustment

The serial CPU block adjustment algorithm of Vlasiator is described as a flowchart in Figure 5. In section a), for each spatial cell, all existing velocity blocks are analyzed. If any of the velocity cells within the block have $f_{s}>f_{\mathrm{th}}$, the block is considered a content block. The loop over the block contents can end early once the first value above the threshold is found. Blocks are thus categorized into two complementary lists (vectors, black lettering) as either being classified as content blocks (green) or non-content blocks (red). Section a) must be completed for all spatial cells before moving on.

Next, in section b), the spatial and velocity halo requirements are fulfilled. For each cell, all its spatial neighbours are identified. For all these spatial neighbour cells, their lists of content cells are looped over, with all contained GIDs added to a required content std::unordered_set (white lettering). For the local cell in question, any local blocks with content are also looped over, and all these blocks and their 26 velocity neighbours are also added to this set. The set contents are indicated as blocks with content (green), blocks with no content (purple), and non-existing blocks (blue, yet to be added).

The third phase, section c), involves pruning non-required blocks from the velocity mesh. The list of local blocks without content is looped over in reverse order (that is, starting from the largest GID i.e.from towards the end of the allocated buffer). Any blocks which are not found in the required content set (dark red) are deleted from the mesh. If the block was the last entry, the mesh size is simply reduced by one block. However, if the block was not the last entry in the mesh, it is replaced with the contents of the last existing block, and then the mesh is reduced in size by one block.

In the final phase, section d), a loop is performed over the required content set. Each GID found in the required content set is compared against the current state of the velocity mesh, appending the velocity mesh with any blocks which do not yet exist (blue). This results in the final state, where the velocity mesh contains all required GID values, each associated with an unique LID forming a consecutive list. This serial approach is straightforward, simple, and results in a relatively low amount of data transfer, but is sequential in nature and does not parallelize well on GPU hardware.

5.2 Parallel GPU block adjustment

The objective of developing the parallel block adjustment algorithm described in this section was to be able to launch GPU kernels which could act in parallel on sections of the same velocity mesh, adding, deleting, adjusting, and moving blocks so that the end result is a valid velocity mesh and so that the kernels do not need to atomically access counters such as taking the next entry from a vector. The algorithm developed for this purpose is displayed in Figure6, and it is run for all spatial cells separately.

In section a), for each spatial cell, all existing velocity blocks (within the localToGlobalMap list or vector, black lettering) are analyzed. The original count $\mathbb{N}_{\mathrm{begin}}$ is stored. A kernel examines all blocks in parallel, and for each block an in-kernel thread-parallel reduction analyzes all contained phase-space cells at once. If any phase-space cell fulfills $f_{s}>f_{\mathrm{th}}$, the block is considered a content block. Blocks are thus categorized into two complementary maps (white lettering) as either content blocks (green) or non-content blocks (red). The maps store the block GID as the key and the LID as the value. After the parallel insertions, the entries from the content map are extracted into a contiguous content list (vector, black lettering). Section a) must be completed for all spatial cells before moving on.

Following the gathering of content lists, section b) ensures all halo requirements are fulfilled. For each cell, all its spatial neighbours are identified. For all these spatial neighbour cells (managed in parallel), their content lists are looped over in parallel, and the GIDs are evaluated against the local content maps. Any neighbour content GIDs found in the no-content map are deleted from there, and instead placed in the content map (purple). If the neighbour GID does not yet exist in either map, it is added to the content map with an invalid LID identifier flagging it as a cell yet to be added (blue). For the local cell, all content blocks and their 26 velocity neighbours are also looped over in parallel, and values are added to the content map and deleted from the no-content map as before. As listed, all of these operations can be performed in parallel. The content map elements are indicated in the figure as blocks with local content (green), blocks with no local content (purple), and non-existing blocks (blue, yet to be added). After all halo requirements have been evaluated, the two resultant maps are evaluated in order to calculate the size of the velocity mesh after all required operations will be complete. The size of the no-content map provides the number of blocks to remove $\mathbb{N}_{\mathrm{remove}}$. The content map is passed to a parallel stream compaction routine, extracting all entries flagged as yet-to-be-added GIDs into a list of blocks to add, and the size of this list provides $\mathbb{N}_{\mathrm{add}}$.

Next, knowing that $\mathbb{N}_{\mathrm{begin}}+\mathbb{N}_{\mathrm{add}}-\mathbb{N}_{\mathrm{$, the generated content and no-content maps can be evaluated in order to prepare lists of GIDs for parallel processing. Using $\mathbb{N}_{\mathrm{after}}$ as an extraction criterion, checking if map entries have a valid value (LID) less than it, we use parallel stream compaction to extract three lists of blocks. Any blocks with content residing beyond the final size of the mesh are designated blocks to rescue. Blocks without content residing beyond the final size of the mesh are to be simply deleted. Blocks without content within the resultant domain of the mesh are to be removed and replaced with either added or rescued blocks.

In the final step, shown in section d) the four constructed lists (blocks to be added, blocks to be deleted, blocks to be replaced, and blocks to be rescued) are used in a parallel operation altering the actual velocity mesh associated with the spatial cell in question. If $\mathbb{N}_{\mathrm{begin}}<\mathbb{N}_{\mathrm{after}}$, the size of the mesh is increased, in order to have sufficient space for all parallel operations. This size adjustment is done on-device if possible in order to avoid page faulting of counters. If the capacity of the map is insufficient, it is performed on host instead. Then, a GPU kernel is launched, where each parallel instance operates on one value or pair of values from the generated lists. Depending on the relative sizes of the four lists, the operation can involve outright deletion of a block, addition and emptying of a new block, replacing an existing block with a new empty block, or replacing an existing block with a block to be rescued. After all these parallel operations have been completed, the contents of the velocity mesh are in their final state. Then, if $\mathbb{N}_{\mathrm{begin}}>\mathbb{N}_{\mathrm{after}}$, the size of the mesh is decreased to match the contents.

6.1 Time step constraints

Although Vlasiator’s semi-Lagrangian solvers do not place a CFL condition as such, a similar limitation is in place due to a limited thickness of the spatial ghost cell layer. Thus, the dynamic simulation time step may be constrained by the largest advection velocity found in the simulation. To this end, a timestep limit reduction kernel is required, which finds the smallest value of $\Delta x/v_{x}$, $\Delta y/v_{y}$, or $\Delta z/v_{z}$ in the simulation, which, with a safety margin, gives the largest allowed spatial advection time step.

6.2 Plasma moments

Reduction of plasma velocity moments [see, e.g., Palmroth etal., 2018] for the simulation domain is performed by gathering device-side pointers to the velocity mesh and block container structures for all local spatial cells, preparing output buffers, and then calling a massively parallel reduction kernel over all local spatial cells. Earlier versions of Vlasiator called a reduction kernel on one cell at a time, but the overhead from repeated kernel calls and back-and-forth transfer of minuscule units of data was prohibitively slow. Kernel merging with no other changes was found to provide a speed improvement of order $10\times$.

6.3 Output data reducers

A single full state of the Vlasiator simulation, as seen in a simulation checkpoint file, can be several terabytes. In order to facilitate data analysis with more reasonable file sizes, output files with reduced data are generated at regular intervals. The most commonly used values stored in these files are values such as electric and magnetic field vectors and the aforementioned plasma moments. However, more detailed data reducers operating over the velocity space of each spatial cell, such as measures of non-Maxwellianity or sub-region analysis (particle precipitation, moments of supra-thermal sub-populations) are also available. Vlasiator implements a flexible ARCH-lambda approach for data reduction, where the calculations to be performed in order to extract the necessary information are prepared as a C++ lambda function, which is passed over to an interface which then either passes it to CPU-based loops, or if GPU hardware is being used, constructs a kernel to which it is passed. The lambda approach as such does not accrue noteworthy overhead, but if the reduction kernel is called for each spatial cell in turn, API overhead as described in Section 6.2 can become prohibitive.

8.1 Simulation setup

The simulation test case consists of a plasma shear flow interface, with two counter-propagating domains of magnetized plasma. The simulation domain is compact and 2-dimensional, spanning 32 cells in the $X$-direction and 28 cells in the $Y$-direction with cell edges of extent 500 km. In the left ($-X$) half of the simulation domain, plasma has a density of $1\,\mathrm{cm}^{-3}$ and a bulk flow velocity of $V=(0,-187.5,0)\,\mathrm{km\,s}^{-1}$. In the right ($+X$) half of the simulation domain, plasma has a density of $2\,\mathrm{cm}^{-3}$ and a bulk flow velocity of $V=(0,+187.5,0)\,\mathrm{km\,s}^{-1}$. The magnetic field is uniform with a value of $B_{z}=12\,\mathrm{nT}$, resulting in an ion gyroperiod of $5.4268\,\mathrm{s}$. The domains are set to have a pressure balance, with a uniform pressure of $22.3\,\mathrm{nPa}$, resulting in temperatures of $384\,\mathrm{MK}$ in the left half and $192\,\mathrm{MK}$ in the right half of the simulation domain. The central transition follows a tanh profile with a transition width of $1000\,\mathrm{km}$. In addition, a small velocity perturbation of amplitude $3.5\,\mathrm{km\,s}^{-1}$ is superposed, with a wavelength spanning the simulation box in the $Y$-direction. In each spatial cell, the velocity space is constructed with maximal extents of $\pm 3000\,\mathrm{km\,s}^{-1}$ in each direction and a velocity resolution of $12.5\,\mathrm{km\,s}^{-1}$ in all of $\Delta v_{x}$, $\Delta v_{y}$, and $\Delta v_{z}$. Figure 7 shows the simulation states at initialization and after $300.2$ s of simulation, showing how a Kelvin-Helmholtz instability [Tarvus etal., 2024] develops causing uneven plasma mixing at the interface. The simulations are run for a total of $700$ s which corresponds to almost 129 gyroperiods, with the KH instability growth saturating after about 50 gyroperiods [as shown in Figure 4 of Tarvus etal., 2024]. This simulation time results in 2337 timesteps when using a value of WID=4, and 2494 timesteps when using a value of WID=8. This difference stems from WID=8 resulting in a stricter CFL translation constraint due to larger velocity space extents, which in turn arise from the large ghost halo domains, as exemplified in Figure 4. We note that this simulation usefully has a stable computational cost (average block counts remain roughly constant) throughout whilst still involving two different plasma domains and showcasing the relevant dynamics.

8.2 Performance comparison on NVIDIA hardware

Both CPU and CUDA versions of Vlasiator were used to run the aforementioned test simulation. For both architectures, variants with both WID=4 and WID=8 were compiled and simulated. The WID=4 versions resulted in approximately 52M blocks in the simulation in total, so roughly 80K blocks per spatial cell. This is greater than usually seen in Vlasiator simulations, but a realistic count in e.g.hot magnetosheath portions of 2D-3V magnetospheric simulations. By comparison, the WID=8 versions resulted in approximately 8.8M blocks in the simulation in total, so roughly 12.7K blocks per spatial cell. The ratio is not quite $8\times$ due to the larger ghost halo domains, as Figure 4 would suggest.

The bar charts in panel a) of Figure 8 show that for CPU variants, most solver portions perform roughly at similar efficiency for both WID=4 and WID=8. A notable exception is block adjustment, where WID=8 takes the lead due to much smaller block lists to manage. The GPU version with WID=4 performs overall somewhat weaker than the CPU versions, due to slightly less efficient block adjustment performance. The total time spent in velocity space updates is larger than CPU variants by almost 80%, which is partly due to the included block adjustment routines, and partly due to insufficient optimization of the small auxiliary kernels of acceleration. With WID=8, however, the GPU version takes a strong lead by utilizing GPU thread parallelism to its fullest. Total propagation time is improved over CPU variants by a factor of three, and both GPU variants show a factor of $>3$ improvement also in spatial translation.

As an additional point of comparison, the simulation was attempted with the GPU version which had the CPU-based serial block adjustment (see Section 5.1) implemented in place of the new parallel GPU block adjustment algorithm (see Section 5.2). This version did however perform the computations on-device in order to reduce page faulting. Nevertheless, this serial approach was much slower than any other version, as shown in panel b) of Figure 8.

These results showcase that with careful design and algorithmic improvements, the capabilities of modern GPU hardware can be harnessed to enhance performance of even complicated simulation simulation codes such as Vlasiator with ever-changing meshes. However, the test case presented here involves complex velocity meshes with large block counts and a relatively small 2-D spatial grid. A similar performance improvement should not be expected in all cases. Preliminary testing with several MPI tasks (not shown) did not yet result in acceptable performance, which we attribute to a lack of GPUdirect RDMA implementation, forcing all MPI communication to go via the host. This hurdle will be tackled in future updates, as efficient multi-node GPU implementations have been shown to be possible [Pekkilä etal., 2022].

8.2.1 Roofline analysis on NVIDIA hardware

In order to gauge the expected future performance improvements which can be coaxed out of Vlasiator’s GPU version, we used NVIDIA’s performance analysis tools to, amongst other things, complete roofline analysis of the Vlasiator spatial translation kernel, as shown in Figure 9. Roofline analysis evaluates the achieved memory transfer rates and the number of arithmetic operations performed by the kernel in relation to the maximum available memory bandwidth and compute capability. The $X$-axis of a roofline plot is the arithmetic intensity, or ratio of arithmetic operations to bytes accessed. The left-hand side of the plot shows a slope, indicating maximum possible performance as limited by available memory bandwidth. On the right-hand side of a roofline plot operations are limited by arithmetic compute capacity. The Vlasiator translation kernel contains both double and single precision floating point operations and has an arithmetic intensity of roughly 3, placing it in the memory-bound regime. The $Y$-axis indicates actual performance in floating point operations per second (FLOP/s) with the Vlasiator translation kernel achieving on the order of 250 GFLOP/s.

We note that this roofline analysis was somewhat artificial due to constraints placed by the NVIDIA CUDA toolchain. Instead of using the usual per-thread per-stream temporary buffers as mentioned in Section 4.3, this test was performed with only a single thread and stream. The temporary buffer of this thread/stream was ensured to be sufficiently large to allow the whole translation task to fit in memory at once. Even so, analysis of this profiling indicated that the kernel experienced stalling during memory reads, explaining why the achieved performance is roughly an order of magnitude below the maximum achievable value. A significant portion of these stalls were long scoreboard stalls when loading VDFs from global accelerator memory. Future optimizations to memory access patterns and pre-caching may improve the performance of this section.

8.2.2 Overall performance on CUDA and NVIDIA hardware

The performance of the acceleration portion of Vlasiator’s solvers on GPU hardware is less straightforward to evaluate, as it consists of several small kernels to be operated in succession. Additionally, the acceleration process updates spatial space for one cell at a time (per thread and stream), so filling the whole device with one kernel invocation is not possible, barring a comprehensive re-write of the whole acceleration framework. This may be a task for a future code update.

In addition to per-kernel performance, total Vlasiator GPU performance was found to be impacted considerably by page faulting and GPU API calls (see Section 2.4) due to several small kernel invocations. This is exemplified in Figure 10, which shows a timeline over one simulation timestep for the GPU WID=8 run as visualized via NVIDIA Nsight Systems, for the previously presented KHI simulation. We utilize the phiprof [von Alfthan, 2023] profiling library, which involves user-defined nested timer regions building up to a profile tree. When compiled on a GPU-compatible platform, phiprof can leverage NVTX (NVIDIA) or ROC-TX (AMD) tracing hooks, allowing the same nested timer ranges to show up in GPU-specific profiling tools.

As an example, the first section of block adjustment (“Com…” within the “re-adjust blocks” timer range) in Figure 10 consists of computing content lists for all spatial cells, a task which has already been optimized to minimize page faults. Before this optimization, page faulting of small counters resulted in spending roughly $3\times$ the amount of time in this range. A future optimization task is to consolidate content list updates into a few large kernel launches which is expected to provide an up to $10\times$ performance improvement, based on the results of Section 6.2.

8.3 Performance on HIP/ROCm and AMD hardware

To evaluate the performance of Vlasiator on AMD hardware, and to validate the usefulness of the unified codebase approach of Vlasiator, we also ran the same Kelvin-Helmholtz simulation on LUMI-G utilizing HIP/ROCm and AMD hardware. To speed up access to unified memory allocations on LUMI-G, we pre-loaded a custom memory manager (https://github.com/sfantao/vlasiator-mempool). Similar to Figure 10, Figure 11 shows a Vlasiator timeline over one simulation time step, showing the different solver sections at work. Comparing the two timelines confirms that overall, the two timelines appear similar, with semi-Lagrangian acceleration taking the longest, and block adjustment and translation taking a quarter to a third of the overall time. On LUMI-G, the translation kernel appears somewhat slower relative to acceleration. This was initially attributed to the low amount of hardware threads resulting in less temporary buffers available for the translation kernel, constraining parallel launch parameters, but increasing the buffer sizes did not appear to alleviate the issue.

This slight relative change in translation versus acceleration is less surprising than the overall poor performance on the tested hardware, as shown in Figure 12. Despite driving a powerful accelerator with several CPU cores and hardware threads, LUMI-G performance was roughly 3.7% of NVIDIA performance with WID=8 and 5.9% with WID=4. A single AMD MI250x GCD is not expected to match the performance of a single NVIDIA A100 card, but they should be able to achieve the same magnitude of performance. Even the optimized straightforward reduction kernels used for velocity moments of the VDFs (see Section 6.2) showed a $50\times$ penalty in performance. This may be partially due to device drivers in need of updating (which is the responsibility of system administrators, not end users) or, for example, slower response time in the unified memory paging subsystem (see Section 4.1). A positive note is that running several MPI tasks was able to benefit from the increased throughput, with the LUMI-G system supporting direct network access from GPU memory (GPUDirect RDMA). The apparent strong scaling, going from 1 MPI task to 8 tasks, is not representative of overall scaling as a single MPI task will be able to skip all network communication entirely. Additionally, the presented simulation case was chosen to have large block counts, which results in heavier-than-usual MPI communication patterns. Multi-node GPU performance analysis will be the topic of a subsequent study.