Barnes–Hut auf Dreiecken

Barnes–Hut ersetzt viele weit entfernte Quellen durch wenige Massenaggregate. Die Idee benötigt keine Würfel: Hier wird sie auf einen binären Baum ähnlicher Dreiecke übertragen. Die Form des Gebiets und die physikalische Massenverteilung bleiben dabei zwei verschiedene Dinge.

1. Warum überhaupt ein Baum?

Die direkte Beschleunigung eines Teilchens benötigt Beiträge aller anderen Teilchen. Für N Ziele und N Quellen wächst die Arbeit quadratisch. Barnes und Hut schlugen 1986 eine hierarchische Gruppierung vor: Nahe Quellen werden detailliert besucht, entfernte Gruppen gemeinsam angenähert [1]. Ihr ursprünglicher räumlicher Baum benutzt kubische Zellen; unser Projekt übernimmt das Gruppierungsprinzip auf der Dreieckshierarchie.

Die geometrischen Dreiecke begrenzen, wo Quellen liegen können. Sie werden für Teilchenkräfte nicht als homogen mit Masse gefüllte Körper interpretiert. Jede Blattquelle behält ihre tatsächliche Position. Für Hydro ist die Blattquelle dagegen eine ausdrücklich gewählte Zellpunktquadratur.

2. Das verwendete Kraftgesetz

\[\boldsymbol a_i=-G\sum_{j\ne i}m_j\frac{\boldsymbol x_i-\boldsymbol x_j}{(|\boldsymbol x_i-\boldsymbol x_j|^2+\epsilon^2)^{3/2}}.\]

G ist die Gravitationskonstante, mⱼ die Quellmasse und ε die gemeinsame Glättungslänge. Die Bewegung erfolgt in einer Ebene, der Kern stammt aber aus dem dreidimensionalen Newton-Potential. Eine echte zweidimensionale Poisson-Gleichung hätte einen logarithmischen Fundamentalkern und wäre ein anderes Modell.

Für große Abstände geht der geglättete Kern gegen die Newton-Kraft. Für r≪ε ist die Kraft proportional zu r statt singulär. ε ist deshalb ein Modellparameter, kein bloßer numerischer Reparaturschalter. Eine Auflösungsstudie mit gleichzeitig veränderter Glättung vergleicht nicht mehr nur denselben räumlichen Diskretisierungsfehler.

\[W=-G\sum_{i<j}\frac{m_i m_j}{\sqrt{|\boldsymbol x_i-\boldsymbol x_j|^2+\epsilon^2}}.\]

Die Paarenergie zählt jedes Paar einmal. Im isolierten kontinuierlichen System sind Energie, Impuls und Drehimpuls erhalten. Daraus folgt noch nicht, dass eine approximative Baumkraft oder ein endlicher Zeitschritt diese Größen exakt erhält.

3. Was ein innerer Knoten speichert

Aus den Kindern L und R berechnet der Aufbau die Gesamtmasse M und den Massenschwerpunkt c. Leere Gruppen werden übersprungen.

\[M_S=M_L+M_R,\qquad \boldsymbol c_S=\frac{M_L\boldsymbol c_L+M_R\boldsymbol c_R}{M_S}.\]

Der geometrische Schwerpunkt des Dreiecks hängt nur vom Bitcode ab. Der Massenschwerpunkt hängt von den Quellen und ihren Massen ab. Er kann nicht aus der Adresse allein gelesen werden.

Warum ist der Monopol sinnvoll? Schreibe xⱼ=cₛ+δⱼ. Durch die Definition des Massenschwerpunkts verschwindet die massengewichtete Summe der Verschiebungen. In der Taylorentwicklung eines glatten Fernfeldkerns entfällt damit der lineare Term:

\[\sum_{j\in S}m_j\boldsymbol\delta_j=0,\qquad \sum_jm_j f(\boldsymbol r-\boldsymbol\delta_j)=M_Sf(\boldsymbol r)+O(\text{zweite Massenmomente}).\]

Für gut getrennte ungeschwächte Quellen liegt der relative Fehler eines einzelnen Gruppenbeitrags typischerweise in der Größenordnung (D/r)². Das ist keine allgemeine Schranke für die Gesamtkraft: Fast entgegengesetzte große Beiträge können sich aufheben, sodass ein kleiner absoluter Fehler relativ sehr groß wird.

4. Öffnen oder zusammenfassen?

Sei Dₛ der Durchmesser des Quelldreiecks und r der Abstand des Ziels zum Massenschwerpunkt. Das klassische Akzeptanzkriterium lautet:

\[D_S^2<\theta^2|\boldsymbol x_i-\boldsymbol c_S|^2,\qquad D_S^2=2^{3-d(S)}.\]

θ ist ein dimensionsloser Öffnungsparameter. Ein kleines θ verlangt eine bessere räumliche Trennung und führt zu mehr besuchten Knoten. Der Test lässt sich quadriert ausführen; dafür ist keine Quadratwurzel nötig. Die Größe D² folgt aus der Tiefe und kann tabelliert werden, der Abstand zum Massenschwerpunkt nicht.

Schematische Dreiecks-Hierarchie und ihre zielabhängige Öffnung: ferne Quellen werden zusammengefasst, nahe Quellen einzeln besucht. — Berechnetes Lehrbeispiel, keine Produktionsmessung. Links die Quellzellen; rechts die für ein markiertes Ziel akzeptierten Gruppen. Grenzen und Massenschwerpunkte stammen aus derselben Bisektionsregel wie im Text.

Bei θ=0 wird kein innerer Knoten akzeptiert. Man erhält die direkte Blattsumme, möglicherweise in anderer Gleitkommareihenfolge. Für Gas ist das weiterhin die Summe der Zellpunktquellen, nicht das exakte Flächenintegral der Dichte.

5. Eine konservative Integer-Alternative

Für Zielblatt und Quellknoten kann man achsenparallele Begrenzungsrechtecke vorbereiten. Der Abstand δ zwischen diesen Rechtecken ist eine untere Schranke des tatsächlichen Ziel–Massenschwerpunkt-Abstands. Berühren oder überlappen sie sich, ist δ=0.

\[\Delta_\alpha=\max(0,b^-_{T,\alpha}-b^+_{S,\alpha},b^-_{S,\alpha}-b^+_{T,\alpha}),\qquad \delta^2=\Delta_x^2+\Delta_y^2.\]

Das Kriterium D²<θ²δ² ist bei gleichem θ geometrisch strenger als das klassische. Auf einem gemeinsamen dyadischen Gitter können die Rechtecklücken und der Vergleich ganzzahlig ausgewertet werden. Wenn θ=p/q exakt den gespeicherten Parameter bezeichnet und 2ʰ die Koordinatenskala ist, lautet der Vergleich:

\[2^{3-d+2h}q^2<p^2\delta_g^2.\]

Exakt ist die geometrische Entscheidung bezüglich der gespeicherten Werte. Die Kraftsumme bleibt Gleitkommaarithmetik. Mehr Öffnungen können kleinere Fehler liefern, kosten aber Zeit. Deshalb müssen Varianten möglichst bei gleicher erzielter Genauigkeit und nicht allein bei gleichem θ verglichen werden.

Die lokale Bitgeometrie-Studie fand keinen Vorteil der unabhängigen Decodierung jedes Knotens beim vollständigen Boxaufbau. Das ist kein Gegenbeweis gegen Bitgeometrie allgemein, sondern ein Befund für diesen konkreten Zugriff und diese Implementierung.

6. Der vollständige Weg zur Kraft

Quellen den Blättern zuordnen
Massentragende Vorfahren aufbauen
Masse und Massenschwerpunkt von unten nach oben aggregieren

für jedes Ziel unabhängig:
    Stapel = [Wurzel]; Beschleunigung = 0
    solange Stapel nicht leer:
        S = entnehme(Stapel)
        falls S leer oder eigenes Blatt: weiter
        falls S innerer Knoten und
              (S enthält Ziel oder Akzeptanztest scheitert):
            lege belegte Kinder auf Stapel
        sonst:
            addiere geglätteten Monopolbeitrag von S

Der nur lesende Baum erlaubt parallele Zielauswertungen. Stark ungleich verteilte Quellen können sehr unterschiedliche Sucharbeit pro Ziel erzeugen; Lastverteilung ist daher ein eigener Messgegenstand. Bei geeigneten Verteilungen ist die Arbeit häufig ungefähr N log N, bei vollständiger Öffnung jedoch wieder quadratisch. Es gibt keine verteilungsunabhängige Geschwindigkeitsgarantie.

Eine einseitige zielweise Baumöffnung erzwingt außerdem nicht exakt gleich große entgegengesetzte Paarbeiträge. Eine gemessene Impuls- oder Energiedrift kann deshalb nicht automatisch allein dem Zeitintegrator zugeschrieben werden.

7. Kraftfehler und Bahnfehler getrennt testen

Die Teilchendynamik verwendet unter anderem Kick–Drift–Kick:

\[\begin{aligned}\boldsymbol v_{n+1/2}&=\boldsymbol v_n+\tfrac{h}{2}\boldsymbol a(\boldsymbol x_n),\\\boldsymbol x_{n+1}&=\boldsymbol x_n+h\boldsymbol v_{n+1/2},\\\boldsymbol v_{n+1}&=\boldsymbol v_{n+1/2}+\tfrac{h}{2}\boldsymbol a(\boldsymbol x_{n+1}).\end{aligned}\]

Mit glatten exakten Hamiltonschen Kräften und festem Schritt ist das Verlet-Verfahren zweiter Ordnung und symplektisch, aber nicht exakt energieerhaltend. Diese Eigenschaften übertragen sich nicht ungeprüft auf jede approximative Baumkraft. Zuerst wird die Kraft gegen die direkte Summe geprüft, danach der Zeitschritt bei festgehaltener Kraftapproximation und schließlich ihre gemeinsame Wirkung.

Quellen und passende Versuche

Barnes & Hut (1986): A hierarchical O(N log N) force-calculation algorithm, Nature 324, 446–449.
Barnes (2012): Gravitational softening as a smoothing operation.
Projektformeln und Einschränkungen: Arbeit 0.3, Kapitel 4 und 7.
G4: direkte Kraftreferenz · Bitgeometrie-Pilotserie · Zeitintegration.