barnes_hut_solver.hpp File Reference

The hierarchical solver: the same physics as direct summation, with most of the interactions replaced by an approximation whose error is bounded and measured. More...

#include <span>
#include <string_view>
#include <vector>
#include "orrery/backend/executor.hpp"
#include "orrery/backend/worker_statistics.hpp"
#include "orrery/core/aligned_allocator.hpp"
#include "orrery/core/softening.hpp"
#include "orrery/core/types.hpp"
#include "orrery/core/vec3_span.hpp"
#include "orrery/solvers/direct_kernel.hpp"
#include "orrery/solvers/force_solver.hpp"
#include "orrery/solvers/interaction_count.hpp"
#include "orrery/solvers/morton.hpp"
#include "orrery/solvers/octree.hpp"

Go to the source code of this file.

Classes
struct	orrery::solvers::EvaluationTimings
	Where the time of one force evaluation went. More...
class	orrery::solvers::BarnesHutSolver
	The Barnes-Hut solver described above. More...

Detailed Description

The hierarchical solver: the same physics as direct summation, with most of the interactions replaced by an approximation whose error is bounded and measured.

Direct summation computes N(N-1) interactions and is exact. This solver groups distant particles into cells and computes one interaction per cell, which turns the cost into something close to N log N and introduces an error controlled by one parameter. That trade is the main algorithmic contribution of this project, and the point of having kept the direct solver is that the error is measured against it rather than argued about.

What actually happens in one evaluation

Four steps, in order, and the first three exist to make the fourth fast.

The particles are sorted along a space-filling curve (solvers/morton.hpp), which puts neighbours in space next to each other in memory. An octree is built over the sorted order, which is cheap precisely because the sort already put every cell's particles in one contiguous range (solvers/octree.hpp). The positions and masses are gathered into that order, so that the traversal reads them sequentially rather than following the permutation on every access. Then every particle walks the tree (solvers/tree_walk.hpp), and the accelerations are written back to the caller's order on the way out.

All four are redone from scratch on every force evaluation. The positions have moved by the time the next one is asked for, and a tree that was updated rather than rebuilt would be a tree whose accuracy depended on how long ago it was built. ADR-0022 records what that costs, which is measured in docs/performance/barnes_hut.md and is a small fraction of the walk.

The accuracy knobs

Two, and they are not equivalent. The opening angle decides how far away a cell has to be before it stands in for its contents, and lowering it buys accuracy by computing more interactions. The quadrupole moment adds the next term of the expansion to every cell that is accepted, and buys accuracy by computing a more expensive interaction rather than more of them. Which is the better way to buy a given error is a question about the machine, and docs/performance/barnes_hut.md answers it by measuring both.

What this solver does not conserve

Momentum, exactly. Direct summation in the form ADR-0015 chose computes each pair from both ends, and the two halves cancel to the last bit because they are the same magnitude computed from the same numbers. A tree evaluation has no such symmetry: particle i may see particle j through a cell while j sees i directly, so the two forces are not equal and opposite, and the total momentum of a configuration drifts.

This is a property of the method rather than a defect of the implementation, it is the reason the direct solver remains the reference, and it is bounded by the same opening angle that bounds everything else. tests/solvers/barnes_hut_test.cpp measures the violation and asserts it stays at the size the approximation implies rather than asserting it is zero, which would be asserting something false.