tree_walk.hpp

Go to the documentation of this file.

#pragma once

 
#include <cstdint>
#include <span>
 
#include "orrery/core/softening.hpp"
#include "orrery/core/types.hpp"
#include "orrery/core/vec3.hpp"
#include "orrery/core/vec3_span.hpp"
#include "orrery/solvers/direct_kernel.hpp"
#include "orrery/solvers/octree.hpp"
 
namespace orrery::solvers {

struct WalkCounts {
    std::uint64_t particle_particle{};
    std::uint64_t particle_cell{};
};

[[nodiscard]] core::Vec3 walk_tree(const Octree& tree, core::Vec3Span<const core::Real> positions,
                                   std::span<const core::Real> masses, core::Index target,
                                   core::Softening softening, AccumulateRange accumulate,
                                   WalkCounts& counts) noexcept;

[[nodiscard]] inline core::Vec3 monopole_acceleration(core::Vec3 offset, core::Real mass,
                                                      core::Softening softening) noexcept {
    // The same expression as one iteration of the direct kernel, with the
    // cell's total mass in place of a particle's and its centre of mass in
    // place of a position. That is what a monopole is, and writing it as
    // anything else would put a second copy of the force law in the project.
    const core::Real factor =
        mass * core::softened_inverse_distance_cubed(core::squared_norm(offset), softening);
 
    return offset * factor;
}

[[nodiscard]] inline core::Vec3 quadrupole_acceleration(core::Vec3 offset, const Quadrupole& moment,
                                                        core::Softening softening) noexcept {
    const core::Real inverse =
        core::softened_inverse_distance(core::squared_norm(offset), softening);
 
    // Formed from the fifth and seventh powers of one reciprocal square root
    // rather than from two further calls, because the divide and square root
    // unit is the bottleneck this project measured in Phase 7 and these are
    // multiplications.
    const core::Real inverse_squared = inverse * inverse;
    const core::Real inverse_fifth = inverse_squared * inverse_squared * inverse;
    const core::Real inverse_seventh = inverse_fifth * inverse_squared;
 
    // Q d, with the tensor's symmetry spent: six stored components describe
    // nine, and the three off-diagonal ones appear twice each.
    const core::Vec3 contracted{
        (moment.xx * offset.x) + (moment.xy * offset.y) + (moment.xz * offset.z),
        (moment.xy * offset.x) + (moment.yy * offset.y) + (moment.yz * offset.z),
        (moment.xz * offset.x) + (moment.yz * offset.y) + (moment.zz * offset.z)};
 
    const core::Real quadratic = core::dot(offset, contracted);
 
    return (contracted * -inverse_fifth) +
           (offset * (static_cast<core::Real>(2.5) * quadratic * inverse_seventh));
}
 
} // namespace orrery::solvers