What is checked, and against what
The code is compared against problems whose answers were known before it was written, and never against its own earlier output except where the claim being made is bitwise reproducibility. Every claim below names the test that makes it, so a reader who does not believe one can run it rather than take it.
That is 311 cases and about eleven seconds. One case skips without an Intel GPU present, and reports that it skipped rather than passing quietly. Four of the presets are checks in their own right rather than conveniences: single-precision runs the same assertions against float with the tolerances that build’s own epsilon implies, sanitise and thread-sanitise run them under the address, undefined-behaviour and thread sanitisers, and sycl adds the GPU cases.
Three references, and which one applies
Analytic results are the strongest statements available, because nothing about them comes from this repository: Kepler’s third law, the virial theorem for a Plummer sphere, the potential energy of a uniform sphere, and the acceleration of a two-body pair at separations where the binary arithmetic is exact.
A compensated reference kernel sums the same softened force law with every term formed in double precision whatever the build’s scalar type, accumulated with Neumaier compensation, so its own summation error is negligible against whatever it is measuring. It is itself checked against the analytic two-body case before it is trusted to judge anything, and that check compares for equality rather than against a tolerance.
Direct summation in double precision is the reference for every approximation: the tree solver at each opening angle, the quadrupole option, the single-precision builds and both GPU solvers. It is not deleted once faster methods exist, and it is what the faster methods have to answer to. Comparing two approximations against each other would say they are wrong in the same way rather than that either is right.
The order of a method, recovered from a measurement
Halving the timestep divides the error of a method of order p by two to the p, so the ratio of the errors at two step sizes recovers p from a measurement. A method whose arithmetic was subtly wrong would still integrate something and would still produce plausible orbits; it would not produce the right power. Measured on a circular orbit, whose exact state after one period is the state it started in.
| Method | Stated order | Measured | Force evaluations |
|---|---|---|---|
| Velocity Verlet | 2 | 1.9998 | 1 |
| Yoshida 4 | 4 | 4.0006 | 3 |
| RK4 | 4 | 4.1670 | 4 |
The same three orders come out of the Python notebooks by an independent route, at 1.9998, 4.0006 and 4.1659, and continuous integration executes those notebooks on every change.
The counterexample
The most informative single comparison the project makes, and the reason it integrates with a second-order method by default rather than a fourth-order one. An eccentric orbit at two hundred steps per orbit, integrated for four hundred orbits, with the energy sampled twenty times per orbit so that the envelope of the oscillation is measured rather than one point on it.
| Method | First twentieth | Last twentieth | Behaviour |
|---|---|---|---|
| Velocity Verlet | 2.6894×10−3 | 2.6894×10−3 | Bounded |
| Yoshida 4 | 2.0170×10−5 | 2.0170×10−5 | Bounded |
| RK4 | 1.92×10−5 | 3.72×10−4 | Growing |
The two symplectic methods return the same error at the end of four hundred orbits as they had at the start, to six digits. A symplectic integrator does not conserve the energy of the system it was given; it conserves the energy of a nearby system exactly, and the difference depends on the timestep rather than on how long the run has been going. RK4 has no such structure to fall back on. Its error is an accumulation, it is nineteen times further out by the end of this run and still moving, and it costs four force evaluations a step against Yoshida’s three to get there.
The same result appears in a column of a CSV file from the command line, over three thousand orbits rather than four hundred.
Velocity Verlet gives 2.685×10−3 in the first twentieth and 2.689×10−3 in the last. RK4 starts twenty times more accurate at 1.30×10−4 and finishes at 2.783×10−3, having just overtaken it. The instrument plays that orbit, and its energy plot is the envelope this table describes.
The two-body problem, and a rotating axis
A two-body orbit has a closed solution, and an eccentric one is sensitive to the exact form of the force law rather than only to its scale. The orbit closes after one period to 8.8×10−6 relative at twenty thousand steps, and to 8.6×10−4 at two thousand. Those two figures are one result rather than two: ten times the step gives a hundred times the error, which is what a second-order method promises and what a mistake in the solver would not deliver.
The precession test is the one worth reading twice. The axis of the eccentric orbit turns by 9.7×10−4 radians per revolution, and a rotating axis is either a force law that is not an exact inverse square, which Bertrand’s theorem makes the only other possibility for a closed orbit, or an artefact of the integrator. The two are told apart by halving the timestep: an artefact falls by four and a wrong exponent does not move. The measured ratio is 3.99. That is velocity Verlet’s second-order phase error and not the physics, and it is the difference between bounding a number and explaining it.
How accurate the kernels are
Against the compensated reference, over a 2048-particle Plummer sphere with a fixed seed. Every available kernel is required to stay inside two hundred times the square root of N times the machine epsilon, and the vector kernel is separately required to be no worse than the scalar one.
| Build | Kernel | Worst | Root mean square |
|---|---|---|---|
| double | scalar | 3.24×10−15 | 8.8×10−16 |
| double | avx2 | 9.5×10−16 | 2.5×10−16 |
| float | scalar | 1.88×10−6 | 4.78×10−7 |
| float | avx2 | 3.60×10−7 | 8.35×10−8 |
Vectorising the sum improves it, by a factor of 3.5 in double precision and 5.7 in single. Keeping one partial sum per lane turns one sum of n terms into four or eight sums of a quarter or an eighth of them, and a shorter sum rounds less. The reassociation a vector kernel is usually apologised for is where some of the error goes rather than where it comes from. No fast-math flag is set anywhere in this project, and the one compiler that defaults to relaxed arithmetic is explicitly put back to strict, because a compiler permitted to reassociate may delete the compensation in the reference kernel and leave the project measuring everything against nothing.
What approximating costs, against what single precision costs
An opening angle of zero reproduces direct summation to 2×10−13, having accepted no cells at all and computed exactly N(N−1) pairs: the same physics reached by a different route, through a different particle order and a different loop. At the default angle the tree’s error is three orders of magnitude larger than the single-precision direct kernel’s own 3.7×10−6 at 65536 particles. Approximating the physics costs far more accuracy than computing it in single precision does, and that is the single most useful accuracy fact the project has measured.
What is not validated
There is no analytic solution to compare a many-body run against. Beyond two bodies the problem is chaotic, so a cluster or a collision is checked through its invariants, its diagnostics and the reversibility of its integrator, and not against a known trajectory. Nothing in this project claims a particular long-duration N-body trajectory is correct in detail, and nothing could.
The tree solver does not conserve momentum, because a tree is not symmetric: particle i may see j through a cell while j sees i directly. The suite measures the size of that violation and requires that closing the opening angle reduces it, which is the true statement, instead of asserting a zero that is not. The renderer’s drawing is not tested either: what is tested is everything under it that is arithmetic, and the drawing itself is checked by looking at it. The GPU evidence comes from one device.
311 cases ·
39 validation
· 36 property
· 3 regression
· 164 unit
Next: every speed, and what it is a fraction
of · play the Kepler run