Orrery
A GPU-accelerated N-body gravitational simulator
Loading...
Searching...
No Matches
diagnostics.hpp File Reference

The conserved quantities, which are how this project knows it is right. More...

Go to the source code of this file.

Classes

struct  orrery::core::Diagnostics
 Every conserved quantity of a configuration, measured at one instant. More...

Functions

Real orrery::core::total_mass (std::span< const Real > masses)
 The sum of the masses.
Vec3 orrery::core::centre_of_mass (Vec3Span< const Real > positions, std::span< const Real > masses)
 The mass-weighted mean position.
Vec3 orrery::core::centre_of_mass_velocity (Vec3Span< const Real > velocities, std::span< const Real > masses)
 The mass-weighted mean velocity, which is the velocity of the centre of mass and is what a sampler subtracts to bring a configuration to rest.
Real orrery::core::kinetic_energy (Vec3Span< const Real > velocities, std::span< const Real > masses)
 T = sum over i of m_i v_i^2 / 2.
Real orrery::core::potential_energy (Vec3Span< const Real > positions, std::span< const Real > masses, Softening softening)
 U = -G sum over pairs of m_i m_j / sqrt(r_ij^2 + eps^2).
Vec3 orrery::core::linear_momentum (Vec3Span< const Real > velocities, std::span< const Real > masses)
 P = sum over i of m_i v_i.
Vec3 orrery::core::angular_momentum (Vec3Span< const Real > positions, Vec3Span< const Real > velocities, std::span< const Real > masses)
 L = sum over i of m_i (r_i cross v_i), about the origin.
Diagnostics orrery::core::measure_diagnostics (const ParticleData &data, Softening softening)
 Every quantity above, measured together.

Detailed Description

The conserved quantities, which are how this project knows it is right.

An N-body system with no external forces conserves its total energy, its linear momentum and its angular momentum exactly. A numerical integration of one does not, and the way it fails is informative: a symplectic scheme keeps the energy inside a bounded envelope for ever while a non-symplectic one of higher order drifts without limit, and that comparison is among the results this project exists to produce. None of it can be shown without these functions, so they are written before the integrators and the solvers that they will judge.

Two properties matter more here than speed.

The potential energy uses the same softened potential as the force kernel (softening.hpp). Softening changes the physics, not just the arithmetic: the system whose energy is conserved is the one with softened masses in it. An unsoftened diagnostic applied to a softened simulation would report a drift that the integrator did not cause, and the conservation result would be measuring the mismatch instead of the method.

The sums are compensated. Total momentum is a sum of terms that cancel to something near zero, and the potential energy of N particles is a sum of N(N-1)/2 terms of one sign; both lose digits under ordinary accumulation, the first to cancellation and the second to a running total that grows far beyond any one term. A diagnostic claiming conservation to round-off has to be more accurate than the thing it measures, so the accumulation is carried with a correction term rather than left to chance.

The potential energy costs N^2 operations, as the direct solver does. It is a diagnostic taken every few hundred steps rather than every step, so this is the right trade: it shares no code with the solver and therefore cannot inherit a bug from it, which is what makes it evidence.

Function Documentation

◆ centre_of_mass()

Vec3 orrery::core::centre_of_mass ( Vec3Span< const Real > positions,
std::span< const Real > masses )
nodiscard

The mass-weighted mean position.

Every sampler in the project recentres on this, so that a configuration's angular momentum is measured about its own centre rather than about whatever offset the sampling happened to produce.

◆ potential_energy()

Real orrery::core::potential_energy ( Vec3Span< const Real > positions,
std::span< const Real > masses,
Softening softening )
nodiscard

U = -G sum over pairs of m_i m_j / sqrt(r_ij^2 + eps^2).

Each pair is counted once. Costs N^2 operations for the reason given at the top of this file.

◆ angular_momentum()

Vec3 orrery::core::angular_momentum ( Vec3Span< const Real > positions,
Vec3Span< const Real > velocities,
std::span< const Real > masses )
nodiscard

L = sum over i of m_i (r_i cross v_i), about the origin.

About the origin rather than about the centre of mass, because that is the quantity an integrator conserves. The two agree for the configurations this project generates, since each is recentred at the origin and left at rest, and where they do not the difference is itself constant.

◆ measure_diagnostics()

Diagnostics orrery::core::measure_diagnostics ( const ParticleData & data,
Softening softening )
nodiscard

Every quantity above, measured together.

The three linear-time quantities share one pass over the particles, which is the shape a simulation wants when it records a diagnostic line; the individual functions above walk the particles once each, which is the shape a test wants when it asks for one number. The potential energy keeps its own loop either way, since it is the only quantity that is not a sum over single particles.