How would you approach a GPU implementation?

[unphased]

I'm delighted to see the parallel branch run at 100fps on my 5800X3D. It's quite fun, and also fun being able to quickly vibe code experiments off of this code base.

GPT-5.2 tells me that it should be pretty practical to do a Barnes Hut GPU implementation. It says:

1. Compute Morton code per particle (based on normalized position).
2. Radix sort particles by Morton code.
3. Build a binary tree over the sorted array using “longest common prefix” logic (Karras 2012 style).
4. Compute internal node AABBs + mass + COM bottom-up.

Why this is nice:

• No dynamic allocation, no atomics-heavy “insert particle into tree” chaos.
• Deterministic structure, fast on GPU.
• Works for 2D or 3D.

Barnes–Hut acceptance test uses each node’s bounding box size and distance, same as quadtree/octree.

On the surface it sounds reasonable. I'll play around with it in my spare time, but I was curious if you had given this any thought and if you think this would be an efficient way to go about it.

Overall (esp in 2026) the quickest way to know if it will be efficient would probably be to just go build it and measure its performance...

[unphased]

After sending it through codex having looked at the code, this further response looks cogent, but will require a few weeks of research for me to make heads or tails of:

• Your suggestion still looks like the right direction for this codebase if the goal is 10–20× more bodies and/or speed.

How it matches (and differs from) what we have now:

• Current code already moved away from “insert one body at a time” toward a data-parallel, allocation-light build by partitioning the bodies array into
  contiguous ranges per node (CPU-side). Your GPU plan is the same philosophy, but taken to the GPU-friendly extreme: Morton keys + radix sort +
  deterministic linear tree (LBVH) instead of recursive partitioning.
• The force evaluation in this repo is “monopole only” (mass + COM), with a Barnes–Hut acceptance test (quad.size^2 < d_sq * theta^2). An LBVH node
  naturally gives you an AABB; you can use s = max_extent (or half-diagonal) and keep essentially the same test: s*s < d_sq * theta^2.

A few practical adjustments I’d recommend when translating your outline to this code:

• Leaves should hold a range of bodies, not singletons (like this repo’s leaf_capacity idea). That keeps traversal depth down and makes the “near-field”
  loop coherent on GPU.
• Morton collisions are real at finite quantization; Karras-style build needs a stable tie-break (often (morton, index)).
• “No atomics” isn’t quite true for bottom-up accumulation unless you do strict level scheduling; the common LBVH pattern uses a tiny atomic counter per
  internal node (“how many children finished”) to know when to accumulate parent mass/COM. That’s still way less chaotic than N inserts with contention.
• Traversal: this repo currently uses a next “rope” pointer for stackless traversal on CPU. On GPU you can either:
  • use a small fixed local stack per thread (often fine), or
  • add a second pass to compute ropes/next pointers for the BVH to keep it stackless.

So yes: Morton → radix sort → Karras LBVH → bottom-up AABB+mass+COM is consistent with the current architecture and likely the cleanest GPU build path,
with the acceptance test mapping directly once you decide what “node size” means from an AABB.