sycl_direct_solver.hpp File Reference

Direct summation on the integrated GPU. More...

#include <memory>
#include <span>
#include <string_view>
#include "orrery/backend/sycl_device.hpp"
#include "orrery/backend/worker_statistics.hpp"
#include "orrery/core/softening.hpp"
#include "orrery/core/types.hpp"
#include "orrery/core/vec3_span.hpp"
#include "orrery/solvers/force_solver.hpp"
#include "orrery/solvers/interaction_count.hpp"

Go to the source code of this file.

Classes
struct	orrery::solvers::SyclEvaluationTimings
	Where the last evaluation spent its time. More...
class	orrery::solvers::SyclDirectSolver
	Direct summation evaluated on a SYCL device. More...

Detailed Description

Direct summation on the integrated GPU.

The same physics as solvers/direct_solver.hpp, over the same softened potential from core/softening.hpp, in the same units. What differs is where the loop runs and how the sources are staged for it. This is a backend behind the existing solver interface rather than a second solver, which is the distinction section 3 of the implementation plan draws and ADR-0026 records the seam for.

The kernel

One work-item per target particle. Each accumulates the acceleration on its own particle from every source, so the write pattern is the same one ADR-0015 chose in Phase 5: read everything, write only yourself, no reduction between work-items and no atomics anywhere.

Sources are staged through shared local memory a tile at a time. Without that every work-item in a group reads the same source position from global memory on the same instruction, which the hardware coalesces but still fetches repeatedly across the N/tile passes. With it, a group cooperatively loads one tile, synchronises, and then reads it from memory local to the Xe-core. That turns the kernel's global memory traffic from O(N^2) into O(N^2 / tile), which matters because section 2 of the plan establishes bandwidth as the binding constraint on this part rather than arithmetic.

Self-interaction

The CPU kernel excludes a particle from its own sum by splitting the range either side of its index, which costs nothing in a loop. That trick does not survive tiling: a tile is a block of sources shared by a whole work-group, and each work-item has a different index to exclude. So the self term is masked instead. Both the mass and the squared separation are selected, the first to zero and the second to one, because zeroing the mass alone would leave an infinite reciprocal distance in an unsoftened run and produce a NaN from inf * 0 rather than the zero contribution intended.

What it will and will not agree with

Not bit-for-bit with the CPU solver, for the reason the AVX2 kernel is not: the sum is reassociated. Here it is reassociated further, since the order is tile by tile rather than index by index, and the device is entitled to contract a multiply and an add into an FMA. Neither is a relaxation of IEEE arithmetic and ADR-0020's rule against fast-math flags is not weakened.

The larger effect is precision, and not in the way integrated GPUs are usually assumed to behave. The target device reports the fp64 aspect, so a double-precision build constructs and runs rather than being declined. What the aspect does not promise is a useful rate, and the gap between reporting double precision and executing it quickly is exactly the sort of claim this project measures rather than assumes: docs/performance/sycl_direct.md reports both configurations and the ratio between them, which is what decides the precision a run should use.

Accuracy is quoted against the double-precision CPU direct solver in either case, since that is the project's reference. tests/solvers/sycl_direct_solver_test.cpp states the tolerance and derives it from the precision the build was configured with.