Bruno Rothschild, Lennart Mannteuffel

Zur Simulation des kontinuierlichen Prozesses der Gravitation und der Bewegung wird die Euler-Methode verwendet. Dies ist die einfachste Methode, Differentialgleichungen zu lösen. Folgende Annahmen sind werden getätigt: \begin{align} \vec{p}(t+\Delta t) &\approx \vec{p}(t) + \vec{v}(t)\cdot \Delta t \\ \vec{v}(t+\Delta t) &\approx \vec{v}(t) + \vec{a}(t)\cdot \Delta t \end{align}

Diese Methode ist auf längere Sicht fehleranfällig, wie in dem Bild zu erkennen. Die Fehler können mit kleineren Zeitschritten minimiert werden.

Sei $d > 0 \in \mathbb{N}$ die betrachtete Anzahl an Dimensionen.
Sei $\vec{p}_i \in \mathbb{R}^{d}$ die Position eines Planeten i, $\vec{v}_i \in \mathbb{R}^{d}$ seine Geschwindigkeit und $m_i \in \mathbb{N}$ seine Masse.

Für die Planeten mit der Masse $m_i$, $m_j$ und den Positionen $\vec{p}_i$, $\vec{p}_j$ und den Geschwindigkeiten $\vec{v}_i$, $\vec{v}_j$ gilt folgendes für die Berechnung der Beschleunigung.
Sei der Differenzvektor zwischen den beiden Positionen der Planeten $\vec{r}_{ij} = \vec{p}_j - \vec{p}_i$. \begin{align} |\vec{F}_{ij}| &= G \cdot \frac{m_i\cdot m_j}{|\vec{r}_{ij}|^{d-1}} \\ \vec{F}_{ij} &= G \cdot \frac{m_i\cdot m_j}{|\vec{r}_{ij}|^{d-1}} \cdot \frac{\vec{r}_{ij}}{|\vec{r}_{ij}|} \\ &= G \cdot m_i\cdot m_j \cdot \frac{\vec{r}_{ij}}{|\vec{r}_{ij}|^d} \end{align} \[ \vec{a}_{ij} = G \cdot m_j \cdot \frac{\vec{r}_{ij}}{|\vec{r}_{ij}|^d} \]

Existieren mehrere andere Planeten, dann ist die vektorielle Beschleunigung des Planeten i gleich der Summe aller einzelnen vektoriellen Beschleunigungen aller anderen Planeten.

\begin{align} \vec{a}_{i,ges} &= G\cdot\sum_j \left(m_j\cdot\frac{\vec{r}_{ij}}{|\vec{r}_{ij}|^d}\right) \\ &= G\cdot\sum_j \left(m_j\cdot\frac{\vec{p}_j-\vec{p}_i}{|\vec{p}_j-\vec{p}_i|^d}\right) \end{align}

Um plötzlichen nahezu unendlichen Geschwindigkeiten bei $|\vec{r}|=0$ und devision by zero Fehlern entgegenzuwirken kann ein minimales $r$ zwischen zwei Planeten gut sein. \[ \vec{a}_{ij} = G \cdot m_j \cdot \frac{\vec{r}_{ij}}{(|\vec{r}_{ij}| + r_{min})^d} \]

Außerdem ist es hilfreich, anzunehmen, dass alle Körper die gleiche Masse haben, da dies die benötigte Speicherkapazität und die Rechenlast reduziert. Mit $m_i = 1$: \[ \vec{a}_{ij} = G \cdot \frac{\vec{r}_{ij}}{|\vec{r}_{ij}|^d} \]

Eine GPU besteht aus vielen (100, 1000, oder mehr) Prozessoren. In einem Vergleich von einem CPU Prozessor und einem GPU Prozessor ist der von der GPU in Leistung und Funktion dem anderen unterlegen. Die große Menge von Prozessoren auf der GPU ermöglicht denoch eine deutliche höhere Leistung, wenn Aufgaben effizient parallelisiert werden können.
Eine typische RTX 3080 besteht zum Beispiel aus 8704 so genannten „shader units“. Der i9-12900K Intel Prozessor besteht im gegensatz aus 16 Kernen.

• nvidia: https://www.techpowerup.com/gpu-specs/geforce-rtx-3080.c3621
• intel: https://ark.intel.com/content/www/de/de/ark/products/134599/intel-core-i912900k-processor-30m-cache-up-to-5-20-ghz.html

Aber es gibt noch einen Haken. Wenn man auf der GPU rechnen möchte, müssen alle dazu benötigten Daten von der CPU auf die GPU geladen werden.

Auf der CPU könnte die Berechnung folgendermaßen ablaufen:

• Für den Planeten i berechnet man die Summe aller Beschleunigungen der anderen Planeten. Diesen Schritt wiederholt man nacheinander für jeden Planeten, also n mal.
• Dabei eintsteht eine rechnerische Komplexität von $O(n^2)$

Mit der Hilfe der GPU kann diese Berechnung auch folgendermaßen parallelisiert werden:

• Anstatt jeden Schritt nacheinander durchzuführen, werden viele Schritte mit der Hilfe der GPU gleichzeitig berechnet.
• Um nicht unnötige Bandbreite der GPU zu blockieren werden zu Beginn alle Positionen, Geschwindigkeiten und Massen in den Speicher der GPU geladen und nur noch dort verändert. Mit Lesezugriffen auf die GPU können Daten abgerufen werden. Alternativ können die Daten auch direkt auf der GPU gerendert werden, was den Datentransfer zwischen GPU und CPU auf ein Minimum beschränkt.

Es gibt verschiedene Sprachen mit denen die GPU programmieren kann. Für Nvidia GPU's gibt es z.B. die hardwarenahe Sprache CUDA C/C++. Auch wenn es mit Einschränkungen verbunden ist, kann diese auch direkt in Python verwendet werden. Die Bibliothek Numba enthält dafür eine Api.
Viel flexibler sind jedoch die GPU-Hersteller unabhängigen Sprachen, wie z.B. HLSL (High Level Shading Language, DirectX, C#-like Syntax) und GLSL (OpenGL Shading Language, OpenGL, C\C++ Syntax). Um GLSL zu benutzen zu können muss auch die OpenGL Graphics Pipline verwendt werden. HLSL hingegen kann man bei der Unity Engine relativ einfach benutzen. Wir haben uns für den Cuda Wrapper Numba und die Unity Engine entschieden.

Als erstes wurden die obigen Grundlagen in einem einfachen Python-Programm mit numpy und matplotlib umgesetzt. Dieses konnte bis zu ca. 200 Körper in einem 2D-Raum simulieren. Da die Berechnung aber erfolgreich war, wurde als nächstes die Beschleunigung der Berechnungen durch Verschiebung auf die Graphikkarte verfolgt.

1000 Objekte, nicht RT

Um von der Grafikkarte zu profitieren zu können, haben wir als nächstes eine komplexere Version in Python erstellt. Dazu wurden die folgenden Bibliotheken verwendet:

• Numba - für Grafikkartenberechnungen
• Pygame - für die Darstellung in einem Fenster
• OpenCv - für die Verarbeitung von Videos
• Pillow - für die Verarbeitung von einzelnen Videos
• Numpy - zum erstellen von Anfangsparametern (z.B. Planetenpositionen)

Dann wurde eine Klassenhirachie erstellt, die es ermöglichen soll die Simulation in Echtzeit auszuführen:

• GameObjekt - Interface für Objekte die mit dem „Game“ interagieren wollen
• Game - organisiert alle GameObjekte und kontrolliert den Main loop (Pygame) und speichert alle Frames in einem Video (OpenCv)
• Universum - kümmert sich um die Berechnung auf der Grafikkarte (Numba), speichert einzelne Bilder (Pillow)
• UniversumGameObjekt - ermöglicht es ein Universum's Objekt im Game zu rendern

Erwartung: Die neue Version der Simulation sollte deutlich performanter laufen als die ohne Grafikkartenbeschläunigung.

Ergebniss: Die Simulation läuft bei 100 Körpern bei ca. 15-20 Fps bei einer Auflösung von 400×400. Die Auflösung hat einen großen Einfluss auf die Performance, da ein Shader verwendet wird, der den Abstand von jedem Planeten zum jeweiligen Pixel berechnet. Dh. es wurde hauptsächlich eine deutliche optische Verbesserung erzielt. Einige Bilder und Videos aus der Simulation:

1000 Objekte, Pixelshader

100 Objekte, experimenteller Pixelshader

100 Objekte, experimenteller Pixelshader

Verbesserungsideen: Vermutlich wird von Pygame das vom Pixelshader erzeugte Bilder ein zweites mal auf dem Prozessor gerendert. Der Pixelshader bietet auch Optimierungspotenzial. Das Programmieren von Shadern ist mit Numba sehr aufwändig.

Die Rechnung wurde in einem Unity Compute Shader implementiert. Damit lässt sich der rechenintensivste Teil des Programmes relativ einfach auf die Graphikkarte verschieben. Im gleichen Compute Shader wird eine Textur erstellt, auf die die Körper dann als einzelne Pixel bzw. den Pixelshader aus den Python-Programmen verwendend, in einer orthographischen Projektion gerendert werden. Dieser ganze Prozess ist sehr ressourcenfreundlich und lässt das Programm bis zu ca. 70000 Körper flüssig in RT simulieren. Eine Drehung der Projektion wurde ebenfalls implementiert, um die Dreidimensionalität sichtbarer zu machen. Als Startbedingungen wurden hier eine Zufallsverteilung der Punkte innerhalb einer Kugel um den Koordinatenursprung und eine laterale Bewegung beider Hälften der Kugel in unterschiedliche Richtungen, um dem System Rotation zu verleihen, verwendet.

70000 Objekte, Pixelshader, nicht RT, Farbe kodiert nach Richtung der Geschwindigkeit

Eine VR-Version von dem obigen Programm wurde zusätzlich erstellt. Dafür wurde die Unity interne Graphikengine verwendet. Die Daten von der GPU-Berechnung werden hier jeden Frame an die CPU gesendet. Dann setzt die CPU die Position von Kugeln auf die berechneten Punkte. Dabei nutzen die Kugeln GPU-Instancing, was die Rendering- und Speicherlast reduziert. Der zeitintensivste Teil dieses Programms ist die Übertragung der Daten auf die CPU, die von der Datenmenge, also von der Anzahl der Körper abhängt. Deswegen kann dieses Programm flüssig nur ca. 30000 Körper simulieren.

Die zweite Version des Unity Projektes nutzt das gleiche Konzept und größtenteils den gleichen Compute-Shader. Der Unterschied liegt im Rendering. In dieser Version wurde ein Shader erstellt, der die Körper im dreidimensionalen Raum mithilfe des eingebauten Renderers rendert. Dies führt zu der Möglichkeit einer beweglichen Kamera. Diese Methode war ursprünglich als effizientere Rendering-Methode für das VR-Projekt gedacht, musste dafür jedoch verworfen werden, da die Kompatibiliät mit dem VR-Renderer innerhalb der Zeit für das Projekt nicht gewährleistet werden konnte.

Python Versionen:

• n-body-simulation-gpu-v1.zip (Numba-Cuda)
• n-body-simulation-gpu-v2.zip (Numba-Cuda)
• n-body-simulation-gpu-v3.zip (Numba-Cuda)

Unity N-Body-Simulation:

• nbody.zip (C#-ComputeShader-HLSL)
• nbodyvr.zip (C#-ComputeShader-HLSL Oculus VR)