Advertisement

Galaxiensimulation mit dem Barnes-Hut Algorithmus

Eine Einführung in das effektive Lösen von N-Körperproblemen

Das N-Körperproblem

Angenommen in einem Raum befinden sich N Körper. Jeder Körper hat ein Potential, das bedeuted er besitzt ein Kraftfeld. Es könnte sich zum Beispiel um eine elektrische Kraft (Coloumbgesetz) oder um eine Gravitationskraft handeln. Die Kräfte der einzelnen Körper überlagern sich zu einem gemeinsamen Kraftfeld. Dieses Feld wirkt wiederrum auf die Körper, die dadurch eine Beschleunigung erfahren (1. Newtonsches Grundgesetz). Wir haben also ein Problem, in dem das Kraftfeld von den Positionen der Körper abhängt und die Position bzw. die Bewegung der Körper wiederum vom Kraftfeld. In aller Regel hat man eine bekannte Anfangsverteilung der Körper mit bekannten Positionen und Geschwindigkeiten und möchte ausgehend davon berechnen, wie sich das System mit der Zeit entwickelt. Typische Anwendungsfälle für derartige Probleme sind zum Beispiel in der Astronomie die Berechnung von Galaxienentwicklung oder in der Biologie das Verhalten von Molekülen.

Die Bewegungsgleichungen sind in aller Regel recht einfach aufgestellt. Wenn man beispielsweise eine Galaxie simulieren möchte, dann lautet die Gleichung um die Kraft auf ein Teilchen der Gesamtmenge zu berechnen:

$$m_i \ddot{\vec{r_i}} = G \sum_{j=1;j \ne i}^N \Big( m_i m_j \frac{ \vec{r_i} - \vec{r_j}} {|\vec{r_i} - \vec{r_j}|^3} \Big)$$

wobei:

Angenommen Fij ist die Kraft zwischen dem Teilchen i und j. Da ein Körper keine Kraft auf sich selbst ausüben kann, ist die Gesamtzahl aller möglichen Kräftepaare im System N*(N-1). Berücksichtigt man weiterhin, das Nach dem Wechselwirkungsgesetz jede Kraft eine gleich große Gegenkraft hat, also Fij = -Fji ist, so reduziert sich die Gesamtanzahl der Kraftberechnungen letztendlich auf:

$${1 \over 2 } N (N-1) \approx O(N^2)$$

Der Aufwand steigt also mit der Ordnung O(N2) an. Wenn sich die Teilchenanzahl verdoppelt, verfierfacht sich der Rechenaufwand, verzehnfacht man die Teichenanzahl, dann verhundertfacht sich der Rechenaufwand. Es ist absehbar, dass diese Art der Lösung des N-Körperproblemes für eine große Teilchenanzahl viel zu aufwendig ist. Möchte man mit sehr großen Teilchenanzahlen rechen, dann benötigt man einen besseren Algorithmus. Daran ändert sich im übrigen auch nichts, wenn man berücksichtigt, dass diese Art von Problem seit einiger Zeit sehr gut auf Spezialhardware (z.B. GPU / CUDA) durch massive Parallelisierung lösen kann. Das Problem verschiebt sich dadurch lediglich um einige Zehntausend Teilchen nach "hinten" und auch der Barnes-Hut Algorithmus kann auf der GPU implementiert werden.

Der Barnes-Hut Algorithmus

Der Barnes-Hut-Algorithmus ist ein 1986 von Josh Barnes und Piet Hut [1] veröffentlichtes Näherungsverfahren, das eine effektive Berechnung der Kräfte in einem N-Körper-Problem ermöglicht. Im Gegensatz zur direkten Aufsummierung der Kräfte, deren Rechenaufwand mit O(N2) ansteigt reduziert sich der Aufwand beim Barnes-Hut-Algorithmus auf O(N log N).

Der Barnes-Hut-Algorithmus verringert die Anzahl zu berechnenden Kräfte durch geeignetes Zusammenfassen von Teilchengruppen zu Pseudoteilchen. Die Grundidee ist dabei, dass die von einer Partikelgruppe auf ein Einzelpartikel ausgeübte Kraft in sehr guter Näherung durch die Wirkung einer einzelnen Masse im Massenschwerpunkt der Partikelgruppe approximiert werden kann. So kann man beispielsweise die Kraft, die von der Andromeda-Galaxie auf die Sonne ausgeübt wird, sehr gut durch eine Punktmasse nähern, die sich im Zentrum der Andromeda-Galaxie befindet und die deren Masse hat. Diese Näherung ist allerdings nur zulässig, wenn der Abstand der Gruppe vom Einzelteilchen groß und der Gruppendurchmesser im Verhältnis zum Abstand klein ist. Das Verhältnis von Gruppenentfernung (d) zu Gruppendurchmesser (r) wird als Multipol-Akzeptanz Kriterium (engl.: MAC; Multipole-Acceptance-Criterion) bezeichnet:

$$\theta = {d \over r }$$

Unterschreitet θ eine bestimmten Schwellwert, so sollte die Näherung nicht mehr angewendet werden um größere Fehler zu vermeiden. Barnes und Hut verwendeten in ihrer Simulation einen Wert von θ=1. Prinzipiell gilt: Je kleiner θ, um so genauer werden die Simulationsergebnisse. Der Barnes-Hut-Algorithmus wendet dieses Prinzip rekursiv an, so dass es auf verschiedenen Längenskalen wirkt (Bild 2). Wie genau das funktioniert wird im nächsten Abschnitt erläutert.

Bild 1: Die von einer Partikelverteilung auf ein Einzelpartikel ausgeübte Kraft kann durch die von einer Punktmasse ausgehenden Kraft angenähert werden, wenn die Entfernung der Partikelgruppe zum Einzelteilchen groß im Verhältnis zum Gruppendurchmesser ist.
Bild 2: Die gravitative Wirkung von Sternhaufen und Stern B auf Stern A kann infolge der großen Entfernung als Punktmasse approximiert werden. Doch auch innerhalb des Sternhaufens kann die Gravitationskraft des Sternhaufens auf Stern B durch eine Punktmasse angenähert werden, da Stern B weit genug vom Sternhaufen entfernt ist.

Die Implementierung

In den folgenden Kapiteln wird die Implementierung des Algorithmus erläutert. Alle Erklärungen und Beispiele beziehen sich auf den zweidimensionalen Fall. Eine Verallgemeinerung auf die dritte Dimension ist unkompliziert, soll hier aber nicht näher erläutert werden.

Die Baumstruktur

Ausgangspunkt ist eine gegebenen Verteilung von Teilchen im zweidimensionalen Raum. Der erste Schritt des Barnes-Hut Algorithmus ist das einsortieren der Teilchen in eine hierarchische Struktur, den Quadtree (bzw. Octtree in 3D). Dafür wird der Raum in vier Quadranten unterteilt. Diese werden im folgenden nach den Himmelsrichtungen NW, NO, SW, SO benannt. Um ein Teilchen in den Baum einzusortieren wird ermittelt, in welchem Quadranten es sich befindet. Wenn sich bereits ein anderes Teilchen in diesem Quadranten befindet, wird der Quadrant erneut unterteilt und beide Teilchen werden in ihre jeweiligen Unterquadranten einsortiert. Dieser Vorgang wird solange wiederholt, bis alle Teilchen einsortiert sind und sich jedes Teilchen in einem eigenen Quadranten befindet. In Bild 3 ist das Ergebnis dieses Vorganges für eine Beispielverteilung dargestellt.

Bild 3: rechts: Räumliche Verteilung der Partikel mit ihren jeweiligen Quadranten; links: Zur Verteilung gehörende Baumstruktur.

Die Erstellung des Quadtree ist ein rekursiver Vorgang, der durch folgenden Pseudocode beschrieben werden kann:

  
Function MainApp::BuildTree
{
  ResetTree

  for all particles
    rootNode->InsertToNode(particle)
  end for
}

Function TreeNode::InsertToNode(newParticle)
{
  if number of particles in this node > 1
  {
    quad = GetQuadrant(newParticle);
    
    if subnode(quad) does not exist
      create subnode(quad)

    subnode(quad)->InsertToNode(newParticle)
  }
  else if number of particles in this node == 1
  {
    quad = GetQuadrant(existingParticle);
    if subnode(quad) does not exist
      create subnode(quad)
    subnode(quad)->InsertToNode(existingParticle)

    quad = GetQuadrant(newParticle);
    if subnode(quad) does not exist
      create subnode(quad)
    subnode(quad)->InsertToNode(newParticle);
  }
  else if node is empty
  {
    store newParticle in node as existingParticle
  }

  Increase number of particles
}

Später in der Simulation muss der Baum für jeden Zeitschritt neu berechnet werden. Die Komplexität der Baumerstellung hängt von der Anordnung der Partikel ab (siehe auch [2]). Der Aufwand um ein Teilchen in den Baum einzufügen ist direkt proportional zur Entfernung des Zielknotens vom Baumursprung. Gibt es viele dicht beieinander liegende Körper dauert das Einsortieren länger, da alle Teilchen in einem kleinen Gebiet konzentriert und der Baum dementsprechend oft unterteilt werden muss.

Berechnen der Masseverteilung

Bis jetzt enthält der Baum nur Informationen über die räumliche Verteilung der Körper. Der nächste Schritt in der Simulation ist die Berechnung der Masseverteilung im Quadtree (bzw. Octtree für 3-D Verteilungen). Es geht dabei speziell um die Gesamtmasse aller in den Knoten enthaltenen Körper, sowie die Berechnung ihrer Massenschwerpunkte. Diese Daten müssen für alle Knoten des Baumes, die nicht leer sind, ermittelt werden. Die Berechnung erfolgt rekursiv durch abscannen der Baumknoten. Folgender Pseudocode beschreibt den prinzipiellen Ablauf:

 
Function TreeNode::ComputeMassDistribution()
{
  if number of particles in this node equals 1
  {
    CenterOfMass = Particle.Position
    Mass = Particle.Mass  
  }
  else
  {
    // Compute the center of mass based on the masses of 
    // all child quadrants and the center of mass as the 
    // center of mass of the child quadrants weightes with 
    // their mass
    for all child quadrants that have particles in them
    {
      Quadrant.ComputeMassDistribution
      Mass += Quadrant.Mass
      CenterOfMass +=  Quadrant.Mass * Quadrant.CenterOfMass
    }
    CenterOfMass /= Mass
  }
}

Der Aufwand für diese Prozedur ist O(N log N) [2].

Die Kraftberechnung

Der letzte aber entscheidene Schritt ist die Berechnung der Kraft auf ein beliebiges Teilchen. Dieser Schritt ist es, der den Barnes-Hut-Algorithmus schneller macht als die direkte Aufsummierung der Kräfte. Das wäre jetzt ein passender Moment sich nochmal die Bilder 1 und 2 ins Gedächtnis zu rufen, denn hier wird nun die Implementierung erläutert. Um die Kraft auf ein Teilchen zu berechnen wird der Quadtree erneut rekursiv durchiteriert. Für jeden Knoten wird das Verhältnis θ nach folgender Formel berechnent:

$$\theta = {d \over r }$$ wobei:

Anschließend wird getestet, ob θ unter oder über einem bestimmten Schwellwert liegt. (Meistens verwendet man eine Schwelle von θ < 1). Liegt θ darunter, so kann die Kraft die von dem Quadranten ausgeht genähert werden, indem der Massenschwerpunkt des Knotens und die Knotenmasse in das Gravitationsgesetz eingesetzt werden. Ist θ größer, so muss der Quadrant erneut unterteilt und alle Unterquadranten nach der selben Methode geprüft werden. Iteriert wird solange es noch Unterquadranten gibt. Im ungünstigsten Falle also bis zu den Endknoten, für die ja Schwerpunkt und Masse identische mit dem einzigen darin enthaltenen Teilchen sind. Das bedeuted gleichzeitig, das man für einen Schwellwert von θ=0 die gleiche Lösung bekommt wie bei der direkten Aufsummierung der Kräfte (nur viel langsamer da der Baum ja trotzdem erzeugt werden muss). Die Wahl von θ beeinflußt die notwendige iterationstiefe und damit die Gesamtanzahl der notwendigen Kraftberechnungen. Animation 1 verdeutlicht diesen Zusammenhang exemplarisch.

Abhängigkeit der Anzahl der Kraftberechnungen von der Wahl von theta
Animation 1: Gezeigt wird eine Verteilung von 5000 Körpern. Es werden nur die Knoten angezeigt, die in die Kraftberechnung auf das Teilchen in der Bildmitte eingehen (Koordinatenursprung). Umso grösser θ ist, je weniger Knoten werden für die Berechnung benötigt, umso grösser ist allerdings auch der Fehler des berechneten Kraftvektors. Wie so oft gilt auch hier: Schelligkeit kontra Genauigkeit!

Abschließend noch der Pseudocode zur Kraftberechnung:

 
Function MainApp::CalcForce
  for all particles
    force = RootNode.CalculateForceFromTree(particle)
  end for
end 

Function force = TreeNode::CalculateForce(targetParticle)
  force = 0

  if number of particle equals 1
    force = Gravitational force between targetParticle and particle
  else
    r = distance from nodes center of mass to targetParticle
    d = height of the node
    if (d/r < θ)
      force = Gravitational force between targetParticle and node 
    else
      for all child nodes n
        force += n.CalculateForce(particle)
      end for
    end if
  end
end

Galaxienkollision als Beispiel

Kombiniert man nun den Barnes-Hut-Algorithmus mit einem geeigneten Integrationsverfahren (z.Bsp. Adams-Bashforth), so kann man mit relativ einfachen mitteln in erträglicher Zeit eine 2D Galaxienkollision errechnen:

zurück      weiter

Das könnte Sie auch interessieren: