The tree the Barnes-Hut solver walks, and the moments it carries. More...

#include <cstdint>
#include <span>
#include <vector>
#include "orrery/backend/executor.hpp"
#include "orrery/core/types.hpp"
#include "orrery/core/vec3.hpp"
#include "orrery/core/vec3_span.hpp"
#include "orrery/solvers/morton.hpp"

Classes
struct	orrery::solvers::Quadrupole
	The traceless second moment of a cell about its centre of mass. More...
struct	orrery::solvers::TreeNode
	One cell of the tree. More...
struct	orrery::solvers::TreeParameters
	The choices that decide what the tree costs and how accurate it is. More...
class	orrery::solvers::Octree
	The octree of one configuration, built from Morton-sorted particles. More...

Functions
TreeParameters	orrery::solvers::corrected_parameters (TreeParameters parameters) noexcept
	parameters with the corrections its own documentation describes: the opening angle brought into [0, 1] and the leaf capacity to at least one.

Detailed Description

The tree the Barnes-Hut solver walks, and the moments it carries.

A cell of the tree stands in for the particles inside it. If the cell is far enough from the particle being accelerated, the sum over its contents is replaced by a single term formed from the total mass and the centre of mass, and the N interactions the cell would have contributed become one. That is the whole of the Barnes-Hut idea; everything in this file is the bookkeeping that makes it exact enough to be worth doing and fast enough to be worth building.

The shape in memory

One array of nodes, in the order a depth-first walk visits them, and no pointers anywhere. Two consequences follow and both are the reason for the choice.

The first child of a node is the node immediately after it, so descending costs an increment and reads a cache line the walk has usually already fetched. The second is the escape index in TreeNode::next: the position of the node a walk should continue at when it does not descend into this one, which is the end of this node's subtree. A walk is then a loop over one index with no stack, no recursion and no parent links, and the two branches it can take are ++node and node = next.

The array is also in an order that suits the machine. The particles were sorted along a space-filling curve before the tree was built (solvers/morton.hpp), so consecutive targets walk nearly the same nodes, and those nodes are near each other in this array because the array is in tree order and the tree is in space order.

Nodes are held as structs, and the particles are not

ADR-0004 stores particle components in separate arrays because the force kernel reads positions and masses and never touches velocities, so an interleaved layout would spend bandwidth on data it does not use. The argument does not carry over to nodes, and this file deliberately does the opposite: a walk arriving at a node reads its centre of mass, its mass, its acceptance radius and one of its two index fields, which is all of it. Held as separate arrays those five reads would touch five cache lines instead of one. The layout follows the access pattern in both cases, which is the decision ADR-0004 was actually making.

The quadrupole moments are the exception, and they are in an array of their own for the same reason. They are read only when they are switched on, and they are as large as the rest of a node put together, so carrying them inside the node would make the common configuration pay for the rare one on every node it touches.

What the tree does not do

It is not a spatial index, it answers no queries and it does not survive a force evaluation. Positions change every step, so the tree is rebuilt every step; ADR-0022 records why that is cheaper than maintaining one.

It also holds no particles. A leaf records a range of the sorted arrays, and the summation over that range is done by the direct kernel of Phase 7, unchanged. The tree's job is to make those ranges short and few.

Function Documentation

◆ corrected_parameters()

TreeParameters orrery::solvers::corrected_parameters ( TreeParameters parameters )

nodiscardnoexcept

parameters with the corrections its own documentation describes: the opening angle brought into [0, 1] and the leaf capacity to at least one.

A free function because two places apply it, the tree when it is built and the solver when it is constructed, and a project that wrote the same clamp twice would eventually clamp differently in the two of them.