Four solvers behind one interface

Direct summation on the CPU, a Barnes-Hut tree on the CPU, and both of them again on an integrated GPU in SYCL. They are four implementations of one interface rather than four programs, and which of them is the right choice is a question about the particle count with a measured answer.

Two qualifications belong beside that table. The crossovers are for a Plummer sphere, and a tree’s cost depends on how clustered the configuration is. And they are for eight threads, since the two algorithms do not parallelise equally well. The tree has to reach half again the size to be worth using on the GPU as on the CPU, which is the expected direction: direct summation is the kernel a GPU is built for, and a tree walk is the kernel it is not.

Five decisions the rest of the code rests on

Particles are stored as separate contiguous arrays rather than as an array of structs. The force kernel reads positions and masses and nothing else. Under an array-of-structs layout every cache line it fetched would also carry velocities and accelerations it never touches, wasting a large share of the bandwidth that already binds. Separate arrays also give contiguous vector loads instead of strided gathers.

Virtual dispatch sits at boundaries and never inside a loop. Solvers and backends are selected at run time, so a benchmark or a test can swap implementations from a flag. The cost is one indirect call per timestep ahead of billions of floating-point operations, and it is unmeasurable. No virtual call appears in any loop over particles.

A GPU implementation is a backend behind the solver interface, not a second copy of the solver. Two divergent implementations of the same physics is the usual way a project of this kind decays, so the device sits behind the interface rather than in front of it, and the suite requires the two to agree rather than to be separately plausible.

Precision is selected at build time. The scalar type is double by default and float under a build option. Templating every solver on it would multiply compile times and complicate the SYCL kernels for no practical gain, because a given run is either accuracy-oriented or throughput-oriented and never both.

The direct solver is the reference and is never deleted. Every approximation, whether an opening angle, a multipole order or a reduced precision, is measured against direct summation in double precision. The validation report is that rule spent.

What one interaction costs, and what a tree charges for skipping it

One pairwise interaction is three subtractions, five operations for the squared separation, one add for the softening, one square root, one division, three multiplies and six operations to accumulate: 20 floating-point operations, counting the square root and the division as one each, which is the convention the N-body literature uses and what makes these figures comparable with published ones. A force evaluation reads 4N values and performs 20N(N−1) operations, which at 8192 particles is an arithmetic intensity of 5119 flop per byte against a ridge point of 3.45. Direct summation is compute-bound by three and a half orders of magnitude.

A tree replaces distant groups of particles with one term each, and the interaction counter says something the timings cannot. At 262144 particles the tree computes 78 times fewer interactions than direct summation and is only 20 times faster, so a tree interaction costs about four times a direct one.

Direct summation computes four pairs per AVX2 register from contiguous arrays with no branch and no address arithmetic. A tree’s cell term is a scalar computation at the end of an unpredictable walk, preceded by the comparison that accepted the cell and by the ones that rejected the cells above it. That factor of four, rather than the asymptotics, is what a faster tree has to attack, and it is the reason the GPU traversal was worth attempting at all.

The rest of a tree evaluation is cheap and measured. The Morton sort, the gather and the build together come to under five per cent at the sizes where a tree is the right choice, and the tree build itself is a fifth of one per cent at 262144 particles. That is the evidence for rebuilding the whole tree every evaluation rather than maintaining one.

Where the device is bound, and where it is not

On the GPU the ratio between the divide ceiling and the multiply-add ceiling is 4.1 rather than 27, so the device’s divide unit is comparatively far stronger and the two processors have to be judged against different ceilings. At 1174 Gflop/s the GPU direct kernel is at 42 per cent of the device’s multiply-add ceiling and 17.5 per cent of its divide ceiling, so unlike the CPU kernel it is bound by neither. What takes the remainder has not been identified, and the two candidates are named rather than guessed between: the two work-group barriers per tile, and the local memory traffic the tiling trades global traffic for. Distinguishing them needs a tile size sweep and sub-group instrumentation that no phase has run.

Saying so is the point. Three smaller items are open and each is recorded where it was found: the direct kernel’s tile size and its two barriers per tile have never been swept, the tree’s leaf capacity has never been swept although its leaves average 9.5 particles against a capacity of 32, and the GPU direct kernel’s reproducible throughput dip at 16384 particles has not been explained.

The decisions above are recorded one at a time in the architecture decision record, which is published with the reference documentation.
Back to the method pages · the instrument