1 Introduction↩︎

In 2026, supercomputers rely heavily on general-purpose Graphics Processing Units (GPUs) [1]. The accelerators favour specific types of calculations: they excel when algorithms exhibit regular, contiguous data access and uniform, independent computations. Most importantly, they favour codes that do not stress the host-GPU interconnect, since it remains slow compared to the GPU on-board memory access, which is typically hitting high-bandwidth memory (HBM). This does not qualitatively change with the advent of superchips with tightly integrated, shared memory. Consequently, certain application types are more popular on these systems than others, and data organisation and access patterns become paramount considerations for efficient GPU utilisation.

One algorithmic class not straightforward to port to GPUs are Lagrangian particle methods. Smoothed Particle Hydrodynamics (SPH) is a prominent example of such a Lagrangian method that is widely used for hydrodynamics simulations in astrophysics since the 1970s [2], [3]. It’s valued for its inherent ability to capture flows with strongly localised features that move in space and time. In the present paper, we focus on the Swift astrophysics and cosmology simulation package [4], which uses SPH to create massive-scale simulations of the Universe [5]–[7] using hundreds of billions of particles. Originally developed as a CPU-only solver, porting Swift to GPUs surfaces many challenges symptomatic of Lagrangian methods in general.

Its core compute kernels describe individual particle updates or particle–particle interactions that can be written in an embarrassingly parallel manner, i.e. fit to a GPU. The natural data model for a particle-based code is the Array-of-Structs (AoS). AoS facilitates particle reordering, sorting, and exchange in a distributed memory environment. Moreover, it maps naturally onto how domain scientists conceptualise data. However, AoS is not the preferred layout from a performance standpoint, as fast kernels tend to favour Structs-of-Arrays (SoA). Beyond that, each SPH compute kernel typically reads and writes only a subset of a particle’s attributes. Looping over particles or particle pairs hence yields scattered data access patterns, while offloading stresses memory bandwidth inappropriately as we ship too much data forth and back [8]. Similar challenges arise in other SPH codes (e.g. [9]–[11]).

Performance-aware developers therefore optimise their data layouts, most commonly by switching to a Structure-of-Arrays (SoA) representation. A complete rewrite into SoA is, however, problematic in SPH: optimising for one compute kernel’s access pattern introduces performance penalties in other program phases. Rather than switching persistently, some authors therefore advocate for temporary, on-the-fly data rearrangement. Such transformations can be applied manually [12] or delegated to a compiler [8]. The latter approach, as well as abstraction layers such as Kokkos (cf. work around [13]), is appealing because data storage decisions need not propagate into the compute kernel code. They keep domain-specific calculations free of layout concerns. The precise combination of data formats and reordering steps depends on the exact algorithmic composition, which varies widely with the physics and numerics implemented. Nevertheless, switching from AoS to SoA is well-understood in principle. The present paper argues that switching representations is, as sole data optimisation, insufficient, as it predominantly focuses on compute performance and as support for gather and scatter operations on the devices improves. Instead, significant attention should be paid to the data organisation on the host such that quick data rearrangements and offloading become feasible.

We propose splitting particles such that each logical particle is represented by multiple physical structs, each holding only a subset of all the attributes. This split AoS (sAoS) representation can often be designed so that all particle attributes required by a given kernel are grouped into one physical struct, while attributes irrelevant to that calculation remain untouched in separate structs. Such a split reduces pressure on the memory interconnect. We further propose that this split can be guided by data-type considerations: it is sensible to subdivide structs so that one AoS holds exclusively double attributes, while another stores only floats or integers. This way, each struct access is homogeneous from a data type point of view which suits vector units.

Our results focus exclusively on the host-sided data reorganisation for various formats, yet suggest that this struct-splitting strategy can accelerate the total GPU-offloading time by 10%–25%. However, the work also reveals open questions: it is neither clear how the splitting should be chosen, as the number of potential sAoS variants grows exponentially with the number of attributes. Nor is it obvious whether knowledge of the exact layout should feed back into the source code. Finally, all optimisations presented here are implemented manually.

The remainder of the paper is organised as follows: We start with a brief sketch of our software architecture and an algorithmic blueprint, highlighting the points where data transformation becomes performance-critical. In Section 3, we present the implementation method to address the data access burden. Performance results in Section 4 demonstrate the potential of the approach, before a brief summary and outlook in Section 5 concludes the discussion.

2 Problem Statement↩︎

2.1 SPH Compute Kernels↩︎

In SPH, local fluid quantities at particle positions are mapped onto the particles’ properties as convolutions with a weighting kernel. In discrete form, these become weighted sums over neighbouring particles [14] and give rise to the algorithmic core of SPH in form of interaction loops, i.e. looping over all neighbouring particles, per reference particle, to collect their weighted contributions. The exact number of neighbours depends on factors such as the choice of the weighting kernel, the requested spatial resolution, and local conditions of the particles. In our case, the support radius \(H\) can vary per particle to accommodate large differences in particle number density [15].

Per time integration step, at least two such interaction loop types are performed: First, the density loop determines the local \(H\) neighbourhood and mass densities per particle. Forces acting on the particles can then be computed in a subsequent force loop. Sphenix [16], the SPH flavour used in this work, employs a third gradient loop between the density and the force loops, and other modern SPH formulations might inject even further ones to increase the method’s accuracy and order of convergence. While all kernels read the particle positions, each interaction loop type reads and writes a different overall set of the total number of particle attributes (Table 1).

The pair-wise interaction kernels team up with linear update kernels which simply traverse a set of particles and update some of their attributes. While the interaction loops yield computationally expensive kernels with a local complexity of \(\mathcal{O}(N)\) over the number of particles \(N\) within a neighbourhood, the linear kernels are cheap and ignored from hereon.

2.2 Data Structures and Tasks in Swift↩︎

Like many SPH solvers [9], [11], Swift sorts particles into cells of an octree: The octree is refined up to a leaf cell size of \(2 H_{\mathrm{max}}\), where each cell hosts a set of particles and \(H_{\mathrm{max}}\) is the largest interaction radius of any of its particles. In the vanilla code version, the particles are encoded as a linear sequence, i.e. the particles are held as AoS. The refinement threshold on the cell size ensures any particle’s neighbours needed during interaction loops are contained either within the cell itself or one of its adjacent cells.

This constrains the workload per interaction loop and also ensures a sufficient compute efficiency: For an arbitrary set of particles, only few will actually contribute towards a force or density, respectively, as only few are within the interaction radius. As particles are stored closely within a cell, the likelihood is high that most particles within a cell actually contribute towards the kernel outcome, while particles that are far away and therefore will not contribute are not considered right from the start.

Another advantage of the decomposition of particles over an octree is the fact that Swift can map the deployment of kernels acting on individual cells onto individual tasks. Swift builds up a large task graph describing the SPH scheme, where each individual task represents one step of the algorithmics over a subset of particles. This leads to high scalability on multicore systems [4].

2.3 GPU-Offloading↩︎

While fine-granular tasks over an octree perform well on CPUs, such tasks induce too little work per kernel to benefit from GPU-offloading. Increasing the cell sizes to contain more particles is not a viable option, as particles only interact with a limited number of neighbours. It would increase the workload but lead to unnecessary computations and thread divergence.

Rather than immediately computing a task as soon as it is scheduled, Swift’s GPU-offloading approach gathers ready tasks into larger bundles of work. The bundles of tasks are then deployed to a GPU en bloc (cf. dynamic task fusion in [17]). The offloading hence decomposes into the following steps:

During the packing process, a task’s relevant particle data is pre-processed (e.g. particle coordinates shifted to facilitate periodic boundary conditions) and moved into a buffer on the host (CPU). As each type of interaction loop (density, gradient, and force) requires only a subset of the particles’ data, only the required quantities are buffered.

Once a bundle reaches a sufficient size, its buffer contents are transferred to the GPU, where the interactions are computed. Subsequently, only the strictly required results are transferred back into host-side buffers. We refer to the whole sequence of data transfer, compute plus transfer back as a launch.

Finally, an unpacking process performs all required reductions and transfers the host-side result buffer into the actual particles.

2.4 Target Hardware↩︎

We perform all experiments below on two architectures: A node comprising two Intel Xeon Gold 6430 CPUs (Intel+A30) and a Grace Hopper superchip. The Grace Hopper supports CPU-GPU coupling through fast NVLINK-C2C memory transfer (900Gb/s) compared to the 64Gb/s of the PCIe 4 connected A30.

2.5 The “packing” Bottleneck and GPU Kernel Realisation↩︎

Both the packing and unpacking steps are executed on the host only, while the launch process comprises both memory transfers and kernel execution. Relative fractions of the time spent in offloading cycles for each interaction loop type (density, gradient, force) in each of the three offloading steps (packing, launch, unpacking) for the Intel+A30 and Grace Hopper show how new generations of accelerators diminish the compute time and make the packing and unpacking steps the primary bottleneck [12] (Figure 1). On Intel+A30, the launch step accounts for \(>97\%\) of the total time, whereas the sum of the fraction of time spent in host-side packing and unpacking buffering operations rises from \(\sim 2\%\) to \(\sim 60-80\%\) on the Grace Hopper.

Additional AoS-to-SoA conversations are out of scope here yet would amplify the effect that the GPU kernels start to outpace their memory transfer and packing counterparts. Swift’s current GPU-offloading is realised through native CUDA, but we observe that OpenMP’s map clause realises similar packing and unpacking. Switching to logically shared memory in return implies that all data transfer is delegated to the hardware via the equivalent of cache line misses. While the hardware might be able to overlap memory transfer and computations more efficiently than Swift with its asynchronous memory streams, it implicitly organises data transfers in memory blocks and hence also moves data which is not required by compute kernels in an AoS formulation. In our setups, this makes it slower than explicit memory transfers.

3 Particle Layout Studies and Optimisation↩︎

Given that the GPU compute kernels require AoS data, in the studies presented, packing means converting AoS data into AoS buffers. Since only a subset of each particle’s data is used during each interaction loop, we hypothesize that an efficient packing benefits from removing un-needed struct members within the target buffer, i.e. the compute kernel benefits from spatial proximity of attributes that are actually used by a kernel and no unncessary data are transferred.

The efficiency of such packing hinges upon the data format of the input data. We therefore experiment with various memory layouts for the host data.

3.1 Particle Memory Layouts↩︎

3.1.0.1 Baseline Particle Memory Layouts

In the AoS version, a single particle’s data is held within a single struct. This mirrors Swift’s default AoS memory arrangement. The struct members are taken directly from the Sphenix SPH implementation of the Swift code and are generally arranged in the same order as they are used throughout the code.

To evaluate the influence of the ordering of the members, we also introduce the random-order version, wherein the ordering of the AoS layout’s struct members have been randomly shuffled.

3.1.0.2 Kernel-Driven Particle Memory Layouts

In an optimised data storage scheme, we split the AoS particle into several structs in anticipation of the kernel access patterns. The struct members of these split Array-of-Structs (sAoS) variants correspond to the data required for the three interaction types. However, since a few variables are required in multiple interaction loops, our code offers three subflavours of sAoS variants. The pack-gradient variant prioritises the gradient``packing operation: All variables required for this operation are placed into its associated struct, disadvantaging the force``packing loop’s access thereof. Conversely, the pack-force variant follows the same principle, but prioritises the force``packing operation instead. Finally, the pack-shared variant places all variables shared between those two operations in a separate “shared” struct.

3.1.0.3 Data Type-Driven Particle Splitting

Modern compute units such as AVX vector registers yield their highest throughput if they operate on variables of homogeneous data types for an extended period. Therefore, we additionally investigate whether they can be improved upon by further splitting each of the sAoS structs according to their data types. These variants are labelled as pack-gradient-type, pack-force-type, and pack-shared-type and split their respective baseline further, such that the first substruct hosts only double precision values for example, the second only single precision, then only integers and so forth.

3.2 Particle Data Accessor Methods↩︎

A simple API hiding the memory layout behind getters and setters allows us to swap particle data memory layout variants at compile time. For the different sAoS and SoA memory layouts, this pushes the logic of how to access the correct sub-struct into the API implementation (Algorithm [algorithm:api]), while the compute kernels can be written storage-agnostic. The APIs are developed against a pointer to a struct, which internally is broken down into scattered substruct accesses for sAoS. Different API implementations can be generated from templates providing a specification of data layout.

struct part_data_arrays {      /* holds all global particle data arrays */
  struct part* _part;          /* main particle struct */
  struct subset1* _subset1; /* First subset of particle data */
  struct subset2* _subset2; /* Second subset of particle data */
  /* etc */
};
struct part {                /* main particle struct */
  size_t _index;             /* this particle's index in particle arrays */
  struct part_data_arrays* _global_data; /* access to global p. arrays */
  /* Other variables ... */
};
__attibute__((always_inline)) inline float get_X(const struct part* restrict p) {
  /* Get variable X, where X is stored in a different struct */
  const struct subset1* restrict _subset1 = p->_global_data->_subset1;
  return _subset1[p->_index].X;
}

Working with logic that maps a (logical) struct pointer onto a meta struct (herein referred to as part-struct), which then in turn delegates to content therein is appealing as it allows for a wide variety of (potentially hierarchical) memory layouts. However, it may not be the most efficient access method and it increases the memory footprint per particle (cmp. additional _index and _global_data attributes).

We implement two more versions to be able to assess any access overhead: The global-var accessor variant accesses particle arrays as a global variable. For this variant, each API call is passed a particle’s index within a global arrays instead of the struct part pointer. In the second variant (explicit-var), we explicitly pass a local variable which holds pointers to global arrays as an argument to each API call instead of accessing them as a global variable.

3.3 Loop Fission Variants↩︎

packing and unpacking are computationally challenging as they lack (significant) arithmetics but move and reorder data that is potentially scattered. We therefore investigate whether these two transformations can be improved by modifying and rearranging their innermost loops. The baseline by-particle creates a single loop where the loop body copies all required quantities of a single particle. The by-struct variant splits the loop \(N\) times, where \(N\) is the number of sAoS substructs whose data is accessed by the conversion. Each loop then only accesses a single struct’s data.

In the by-element variant, every single particle quantity being accessed is accessed in a separate loop. Finally, we also consider data type-driven loop fission. With the by-type version, loops are split such that each loop only accesses a single data type. The by-struct-and-type combines this approach with the memory layout driven by-struct approach and splits the by-struct loops further by the accessed variable’s data types.

3.4 Swiftgpupacksim Mini-App↩︎

To zoom in on the packing and unpacking operations, we introduce a purpose-built mini-app “Swiftgpupacksim”¹. It accepts a trace of GPU-offloading actions from a “real” Swift simulation and replays this recorded sequence of packing and unpacking operations. Each thread’s workload is logged separately over 5 full simulation steps, recording realistic data access patterns of \(\mathcal{O}(10^5 - 10^6)\) operations per thread and per simulation step.

The underlying GPU access patterns stem from two different examples. “Gresho256” consists of a Gresho-Chan vortex [18], where a cubic domain with periodic boundaries is discretised using \(256^3 \approx 16.8 \times 10^6\) particles (\(\approx 2.5\)GiB particle data) which are approximately uniformly spaced in a glass-like distribution. The minimal, maximal, mean, and median number of particles per cell (and hence per packing or unpacking operation) are 55, 73, 64, and 64, respectively.

The second used experiment, “Eagle25”, tests the behaviour and performance on highly non-uniform initial conditions extracted from the EAGLE cosmological simulation campaign [6] at late cosmic times, where cosmological structures and voids are well-formed, leading to large differences in particle number densities contained within the simulated volume. It consists of \(\approx 50.7 \times 10^6\) particles (\(\approx 7.74\)GiB of particle data). Here, the minimal, maximal, mean, and median number of particles per cell are 5, 125’113, 184, and 109, respectively.

4 Results↩︎

4.0.0.1 Memory Layout Impact on Overall Simulation

We first assess the total time required to complete all packing and unpacking operations (Figure 2). A smaller runtime indicates a faster conversion kernel, i.e. a better outcome. The
random-order memory layout consistently demands more time (\(\approx 10\%\) on
Intel+A30 and \(\approx 20\%\) on Grace Hopper) than its well-ordered AoS counterpart. With the exception of the Gresho256 experiment on Grace Hopper, all sAoS variants outperform the AoS variant by up to \(\approx 10-20\%\) on Intel+A30 and \(\approx 20-40\%\) on Grace Hopper. Further splitting the sAoS variants by data type generally yields additional, modest runtime improvements of \(\approx 5\%\). These gains are most pronounced when the particle-per-cell count grows very high.

It is intuitively clear that grouping particles together is advantageous, as it allows conversion kernels to transfer contiguous attributes en bloc into the GPU data structure—likely within a single cache line. Splitting data into separate structs (sAoS) yields further speedup: even if substructs are not stored contiguously, the probability that unnecessary struct fields are loaded into cache during a conversion remains low.

Unfortunately, codes such as SPH employ many different compute kernels, so it is not clear from the outset which kernel should dominate the data layout.

4.0.0.2 Impact on individual compute kernels

We next break down the measurements into the impact of data layouts on individual compute kernels (Figure 3) and evaluate whether the access realisation affects conversion kernel performance.

We observe that particular particle sublayouts have a pronounced effect on different kernels: the total impact of a given layout optimisation is therefore a weighted combination of the runtime improvement due to sAoS substructuring and that kernel’s contribution to total simulation time.

The explicit-var and global-var variants show very little difference from each other, while part-struct access is consistently the slowest for packing operations—yet delivers the fastest gradient and force``unpacking. To explain this behaviour, we examined the optimisation reports generated by the compiler (Intel OneApi 2025.2.1). For the part-struct accessor, the compiler is unable to auto-vectorise the innermost loop: the indirection introduced by sAoS obscures any contiguous memory access, even when the data is physically stored as a linearised (sub-)AoS stream. Counterintuitively, this lack of vectorisation sometimes proves beneficial. The outliers in the gradient/unpacking operation (pack-gradient, pack-shared, pack-force-type, and pack-shared-type) are auto-vectorised, yet vectorisation here hurts performance, likely due to excessive scatters and gathers and a reduced core frequency.

4.0.0.3 Loop Fission Variants

Loop fission yields no improvement on Intel+A30; consulting the optimisation reports, the compiler tends to fuse the loops back together. On Grace Hopper, however, splitting loops by-type reduces runtime by \(\approx10\%\) for nearly all memory layouts relative to unsplit loops. The architecture benefits from homogeneous data-type access patterns (Figure 4).

5 Summary and outlook↩︎

If a code rewrites a data structure prior to invoking a compute kernel on a GPU, managing the conversion overhead is key. It is natural to convert only those attributes subject to reads and writes by that particular kernel. In related work, we refer to the subset of affected particle attributes as a data view [19]. The present work suggests that identifying a subset of attributes immediately before an on-the-fly conversion is too late. Instead, we decompose struct data a priori, so that a conversion operates on only a subset of attributes from the outset.

Conversely, our data suggest that GPUs’ compute capabilities improve at a rate such that the textbook wisdom that SoA plays the dominating role is no longer universally true. Indeed, many compute kernels on GPUs are now that fast that data transfer dominates the runtime. This insight challenges traditional AoS-to-SoA conversion techniques [8], but at the same time implies that a proper design of sAoS is important. While the sAoS technique proves valuable in a robust and consistent manner, fundamental challenges remain.

First, our implementation of the data structure variants is manual. This is acceptable for a case study and appropriate for a codebase that solves a single, well-defined problem. For a general-purpose code such as Swift, however, where users are expected to add further data fields to a struct and to inject new computational kernels (e.g. to implement additional physics), manually maintaining data conversions is problematic: It increases code complexity, is error-prone, and complicates the introduction of new features. Therefore, an open question remains of whether the decomposition of AoS into sAoS can be hidden behind an abstraction layer such as Kokkos [20], injected via compiler directives analogous to attributes or pragmas [21], or fully delegated to the compiler.

Second, our work does not investigate whether robust heuristics exist or might exist for choosing the sAoS decomposition. The minimal reasonable granularity appears to be determined by the view, i.e. the attributes’ read and written by a single compute kernel, though a finer subdivision proves beneficial in some cases if it induces conversion loops operating over homogeneous data types. Conversely, fusing the attribute sets of multiple kernels into one sAoS substruct may sometimes be advantageous—for instance, when two kernels inducing different views execute in immediate succession on the same accelerator. Robust heuristics are needed to guide effective performance engineering in this space. Data replication may also be worth considering: if all substructs within an sAoS decomposition require a particle’s position, it may be beneficial to store that information redundantly so that two conversions can proceed independently and without scattered memory access, even though both require the same field.

Finally, our work does not fully resolve the interplay between sAoS decomposition and vectorisation. We expect that loop fission and loop reordering can become powerful tools in combination with sAoS decomposition—particularly in cases where such restructuring allows the compiler to auto-vectorise more aggressively and to overlap the conversion of partial particle data with compute kernel execution. In those scenarios, a physical decomposition of a logical struct may yield significant speedups.

References↩︎