Guide

This guide walks through a complete build-and-search cycle with PhysicalTrees.jl. By the end of it you should be comfortable building a tree (single- or multi-worker), performing every flavour of neighbor search, and tearing the tree down again.

1. Constructing a tree

PhysicalTrees.jl accepts three kinds of input:

A plain Vector{<:PVector} (positions only, no particle metadata).
A StructArray{<:AbstractParticle} (e.g. Star(uAstro), Massless).
A unitless StructArray{<:AbstractParticle} (e.g. Star()).

In all three cases the entry point is octree.

From raw positions

using PhysicalParticles, PhysicalTrees, Unitful, UnitfulAstro

pos = [PVector( 1.0,  1.0,  1.0, u"kpc"),
       PVector(-1.0, -1.0, -1.0, u"kpc"),
       PVector( 1.0,  0.0, -1.0, u"kpc"),
       PVector(-1.0,  0.0,  1.0, u"kpc"),
       PVector( 0.0,  0.0, -1.0, u"kpc"),
       PVector(-1.0,  0.0,  0.0, u"kpc")]
tree = octree(pos)


Tree defined on worker 1:
    Distributed on workers: [1]
                     units: Unitful.FreeUnits{N, D, nothing} where {N, D}[kpc, Gyr, A, K, cd, M⊙, mol]
        Number of topnodes: 953
       Number of topleaves: 834
             Domain factor: 1.0475284715284717e6 kpc^-1
         Domain start list: [1]
           Domain end list: [834]
                 Load list: [6]
                 Work list: [0.0]
      Number of tree nodes: Any[953]
    Extent: xMin = -1.0 kpc, xMax = 1.0 kpc, yMin = -1.0 kpc, yMax = 1.0 kpc, zMin = -1.0 kpc, zMax = 1.0 kpc, SideLength = 2.002 kpc, Center = PVector(0.0 kpc, 0.0 kpc, 0.0 kpc)

  -------------------------- Data info --------------------------
                         total: 6
                          cuts: Any[6]
              last communicate: Dict(1 => 0)

From a StructArray of particles

If you have mass, smoothing length, ID, etc. on each particle, build the particle array first and pass it to octree:

particles = StructArray(Star(uAstro) for _ in 1:6)
assign_particles(particles, :Pos, pos)
tree = octree(particles)


Tree defined on worker 1:
    Distributed on workers: [1]
                     units: Unitful.FreeUnits{N, D, nothing} where {N, D}[kpc, Gyr, A, K, cd, M⊙, mol]
        Number of topnodes: 953
       Number of topleaves: 834
             Domain factor: 1.0475284715284717e6 kpc^-1
         Domain start list: [1]
           Domain end list: [834]
                 Load list: [6]
                 Work list: [0.0]
      Number of tree nodes: Any[953]
    Extent: xMin = -1.0 kpc, xMax = 1.0 kpc, yMin = -1.0 kpc, yMax = 1.0 kpc, zMin = -1.0 kpc, zMax = 1.0 kpc, SideLength = 2.002 kpc, Center = PVector(0.0 kpc, 0.0 kpc, 0.0 kpc)

  -------------------------- Data info --------------------------
                         total: 6
                          cuts: Any[6]
              last communicate: Dict(1 => 0)

From unitless data

using PhysicalParticles
pos_unitless = PVector.([1.0, -1.0, 1.0, -1.0, 0.0, -1.0],
                         [1.0, -1.0, 0.0,  0.0, 0.0,  0.0],
                         [1.0, -1.0, -1.0, 1.0, -1.0, 0.0])
unitless_particles = StructArray(Star() for _ in 1:6)
assign_particles(unitless_particles, :Pos, pos_unitless)
tree_unitless = octree(unitless_particles)


Tree defined on worker 1:
    Distributed on workers: [1]
                     units: nothing
        Number of topnodes: 953
       Number of topleaves: 834
             Domain factor: 1.0475284715284717e6
         Domain start list: [1]
           Domain end list: [834]
                 Load list: [6]
                 Work list: [0.0]
      Number of tree nodes: Any[953]
    Extent: xMin = -1.0, xMax = 1.0, yMin = -1.0, yMax = 1.0, zMin = -1.0, zMax = 1.0, SideLength = 2.002, Center = PVector{Float64}(0.0, 0.0, 0.0)

  -------------------------- Data info --------------------------
                         total: 6
                          cuts: Any[6]
              last communicate: Dict(1 => 0)

2. Distributed construction

The same call works across workers; pass pids = workers() (or any subset thereof) and the build pipeline is split:

using Distributed
addprocs(2)
@everywhere using PhysicalParticles, PhysicalTrees, Unitful, UnitfulAstro

# Re-attach the data on the master so the example block is self-contained
pos = [PVector( 1.0,  1.0,  1.0, u"kpc"),
       PVector(-1.0, -1.0, -1.0, u"kpc"),
       PVector( 1.0,  0.0, -1.0, u"kpc"),
       PVector(-1.0,  0.0,  1.0, u"kpc"),
       PVector( 0.0,  0.0, -1.0, u"kpc"),
       PVector(-1.0,  0.0,  0.0, u"kpc")]
particles = StructArray(Star(uAstro) for _ in 1:6)
assign_particles(particles, :Pos, pos)
tree = octree(particles, pids = workers())
unregister(tree)

Do not run distributed `octree` inside `@example` / `@repl` blocks

Distributed.remotecall_eval plus worker precompilation inside a Documenter-evaluated block can hang the entire docs build. The snippet above is shown as plain code (not evaluated) for that reason. Run it in a real Julia session yourself to see the output.

Under the hood, octree calls init_octree on every worker, then split_domain, then build, then update. All four are exposed for advanced users who want fine-grained control — see the octree page for details.

3. Searching for neighbors

Once you have a tree, the search family is:

Function	Use case
`ngb_treefind_variable`	Single center, fixed cube radius.
`ngb_treefind_pairs`	Single center, fixed sphere radius with per-particle smoothing.
`ngb_treefind_all`	Bulk: neighbor list for every local particle (SPH density loop).
`ngb_treefind_pairs_all`	Bulk version with per-particle `hsml`.
`search_inbox_local` / `search_inradius_local`	Distributed variants that also return a per-worker `ExportFlag`.
`search_all_parallel`	Cube search that automatically dispatches the export to remote workers.

# Cube search: all particles within ±1 kpc of the origin
ngb = ngb_treefind_variable(PVector(0.0, 0.0, 0.0, u"kpc"),
                            1.0u"kpc", tree)
length(ngb)

# Sphere search via the pairs variant
ngb_sphere = ngb_treefind_pairs(PVector(0.0, 0.0, 0.0, u"kpc"),
                                1.0u"kpc", tree)
length(ngb_sphere)

# Bulk neighbor list (used by the SPH density pass)
all_ngb = ngb_treefind_all(tree, 1.0u"kpc")
length(all_ngb), sum(length, all_ngb)

(6, 10)

4. Periodic boundary conditions

For a periodic box of side L, build the tree with OctreeConfig and use the minimum-image search routines:

cfg = OctreeConfig(6; Periodic = true, BoxSize = 2.0)
tree_pbc = octree(particles, pids = [1], config = cfg)
# The particle at (-1, -1, -1) now has the particle at (1, 1, 1) as
# its periodic neighbor when searching near (1.95, 1.95, 1.95):
ngb_pbc = ngb_treefind_variable(PVector(1.95, 1.95, 1.95, u"kpc"),
                                0.1u"kpc", tree_pbc)
length(ngb_pbc)  # contains the (1, 1, 1) particle via the periodic image
unregister(tree_pbc)

5. Moving particles and re-using the tree

The tree structure is built from particle positions at construction time. After a few SPH / gravity steps, particle positions have drifted; the tree's SideLength may no longer contain them. Call update_node_len to enlarge nodes as needed (equivalent to Gadget-2's force_update_node_len_local):

# Move every particle to (100, 100, 100) kpc
new_pos = PVector(100.0, 100.0, 100.0, u"kpc")
for i in eachindex(tree.data)
    tree.data[i] = setproperties!!(tree.data[i], Pos = new_pos)
end

# Grow the tree to contain the new positions
PhysicalTrees.update_node_len(tree)

For a full rebuild (e.g. after a large re-grid), use rebuild:

rebuild(tree)


Tree defined on worker 1:
    Distributed on workers: [1]
                     units: Unitful.FreeUnits{N, D, nothing} where {N, D}[kpc, Gyr, A, K, cd, M⊙, mol]
        Number of topnodes: 953
       Number of topleaves: 834
             Domain factor: 1.0475284715284717e6 kpc^-1
         Domain start list: [1]
           Domain end list: [834]
                 Load list: [6]
                 Work list: [0.0]
      Number of tree nodes: Any[953]
    Extent: xMin = -1.0 kpc, xMax = 1.0 kpc, yMin = -1.0 kpc, yMax = 1.0 kpc, zMin = -1.0 kpc, zMax = 1.0 kpc, SideLength = 2.002 kpc, Center = PVector(0.0 kpc, 0.0 kpc, 0.0 kpc)

  -------------------------- Data info --------------------------
                         total: 6
                          cuts: Any[6]
              last communicate: Dict(1 => 0)

6. Cleaning up

When the tree is no longer needed, free the per-worker PhysicalTrees.registry entries with unregister:

unregister(tree)

If you want to wipe all trees (e.g. before constructing a fresh one with the same id), call clear(workers()) instead.