Kolumne, UTC (aktualisiert UTC)

Newton wusste, dass er falsch lag

Schneller rechnen, Teil 1: Eine Wolke von Punktmassen

Dies ist der erste Text in einer geplanten Reihe zu effizienterem, schnellerem wissenschaftlichem Rechnen. Darin stelle ich eine einfache physikalische Berechnung vor, welche sich als wiederkehrendes Beispiel wie ein roter Faden durch weitere Texte dieser Reihe ziehen wird. Im Rahmen dieses Artikels wird die umzusetzende Modellphysik analysiert, auf dieser Analyse aufbauend eine erste Struktur für eine spätere Simulation vorgeschlagen und diese abschließend optimiert.

In: Wissenschaftliches Rechnen, Berechnung, Simulation, Physik, Optimierung
 

Schneller rechnen

Sie haben ein Programm, eine Simulation oder ein Modell geschrieben. Die Berechnungen sind nicht trivial und überaus zeitaufwendig. Für die Entwicklung blieb nur begrenzt Zeit - Ihre Wahl fiel nicht zuletzt deswegen auf eine Hochsprache wie Python, R, Matlab, Julia, IDL, JavaScript oder Ruby beziehungsweise auf eine Sprache aus Microsofts .Net-Familie. Ihr Produkt funktioniert, aber es ist langsam. Jede scheinbar noch so kleine Erweiterung lässt es noch viel langsamer werden. Sie warten Stunden, vielleicht sogar Tage, auf das Ergebnis eines Durchlaufs und haben sich dafür schon lange den schnellsten Computer in Ihrer Reichweite gesucht. Wirklich akzeptabel für Ihre sonstigen Arbeitsabläufe wären eigentlich nur wenige Minuten. Das Ausführen von Berechnungen auf einem Rechnerverbund beispielsweise in einem Rechenzentrum kommt auf Grund von Datenschutz, Geheimhaltung oder schwer bis gar nicht durchführbarer Parallelisierung in der gewählten Hochsprache nicht in Frage.

Vor diesem oder ähnlichen Problemen stehen heutzutage vermutlich nicht wenige Menschen. Ob im wissenschaftlichen Betrieb oder in Unternehmen - zeitliche Einschränkungen und mangelhafte Bildungsangebote für so genannte "Nicht-Informatiker" führen immer wieder zu solchen Situationen. Programmieren an sich ist dabei eigentlich nicht (mehr) das Hindernis. Vielmehr ist es ein Mangel in der Kommunikation von Wissen um die unzähligen Möglichkeiten, die es heute in diesem Bereich gibt und wie man diese punktuell, sinnvoll und effizient einsetzen kann. An dieser Stelle möchte ich mit einer kleinen Reihe an Texten ansetzen. Sie ist für Leute gedacht, für die Computer nicht den Mittelpunkt ihres Lebens darstellen sondern ein Werkzeug für die Lösung ihrer eigentlichen Probleme sind. Ich möchte Menschen ansprechen, die sich - genau wie ich - ausdrücklich nicht als Programmierer oder Entwickler sondern als Nutzer verstehen, die sich mit den Mitteln moderner Informationstechnologie in ihrem eigentlichen Fachgebiet ausdrücken müssen.

Am Anfang war die Hochsprache

Eine Hochsprache zur Entwicklung von komplizierten Berechnungen zu verwenden ist im Prinzip richtig. Vor allem in der Anfangsphase erlaubt dies schnelles und unkompliziertes Vorankommen. An der relativ geringen Ausführungsgeschwindigkeit ist bis hier her auch nichts auszusetzen - garantieren doch die meisten Hochsprachen im Gegenzug, auf alle Klippen im Unterbau der Berechnung zu achten und diese vor dem Programmierer, der hier eigentlich ein Benutzer kommandoorientierter Software ist, zu verbergen. Gedanken über Angelegenheiten wie Speicherverwaltung entfallen und schaffen Zeit für die eigentlichen Aufgaben, die es zu bewerkstelligen gilt. Diese Strategie wird erst dann wirklich problematisch, wenn an irgendeiner Stelle ein nicht mehr zu vernachlässigender Engpass entsteht, der eine Berechnung weit über das jeweils akzeptables Maß hinaus verlangsamt. Viele moderne Fachbücher, welche Hochsprachen behandeln, raten dieser Stelle neben einigen die jeweilige Hochsprache betreffenden Maßnahmen lapidar dazu, den fraglichen Engpass "einfach" in der Sprache C zu re-implementieren, ihn dort gegebenenfalls zu parallelisieren und ihn so zu entfernen, ohne auf die notwendigen Details einzugehen oder auf geeignete weiterführende Literatur zu verweisen.

In der Realität gestaltet sich dieser Schritt leider selten als einfach. Selbst wenn man über die nötigen Kenntnisse in C verfügt und weiß, wie man ein in C geschriebenes Stück Programm in die eigentliche übergeordnete in einer Hochsprache verfasste Berechnung einbindet, erfordert dieser Schritt in der Regel trotzdem nicht zu unterschätzende Umbaumaßnahmen im gesamten Projekt. Um nicht später in diese Falle zu laufen sollte man bei Zeiten daran denken, dass die Berechnung eines Tages auf die eine oder andere Art beschleunigt werden muss, und dies auch von vornherein bei der Strukturierung der Berechnung einplanen. Auch wenn dies jetzt mehr nach einem mathematischen beziehungsweise rein technischem und damit sekundären Problem klingt, so ist eine entsprechende gute Strukturierung untrennbar mit dem Inhalt der eigentlichen Berechnung verbunden und lässt sich nur mit fundierten Kenntnissen in der jeweiligen Disziplin, der die Berechnung zuzuordnen ist, beispielsweise Physik oder Chemie, bewältigen. Eine entsprechende Strukturierung, welche später eine Beschleunigung erlaubt, sollte also schon bei der Analyse der eigentlichen fachlichen Fragestellung ansetzen.

Eine Wolke aus Punktmassen

Um mich in neue Sprachen und Techniken einzuarbeiten verwende ich schon seit Jahren unter anderem ein relativ einfaches Beispiel: Eine Wolke aus unzähligen Punktmassen, quasi Partikeln, die sich frei durch einen (nahezu) unendlichen Raum, meinem Modelluniversum, bewegen können und sich dabei gegenseitig anziehen. Als Startwert erzeuge ich mir zwei Gebilde aus in etwa gleich vielen Punktmassen, die Spiralgalaxien nicht unähnlich sehen, und lasse diese ineinander fliegen, woraufhin meine zwei Gebilde auf recht beeindruckende, charakteristische Art zerstieben. Um dies zu realisieren, berechne ich in konstanten Zeitschritten die Anziehungskräfte zwischen allen Punktmassen und bewege diese entsprechend der Summe der auf sie wirkenden Kräfte einen Schritt im Raum weiter.

Video 1: Beispielrechnung mit 60000 Punktmassen über 12000 Zeitschritte. Die Berechnung erfolgte auf einem Intel Core2 Quad Q9550 bei 2,83 GHz und benötigte rund 8 Stunden.

Die Grundlage für diese Berechnungen bildet das newtonsche Gravitationsgesetz:

(1)  F = G \frac{m_1 m_2}{r^2}

Dabei ist F der Betrag der Kraft in Newton, die auf beide Punktmassen wirkt, m_1 und m_2 beschreiben die Masse der beiden Punktmassen in Kilogramm, r stellt den Betrag des Abstandes der beiden Punktmassen dar und G die Gravitationskonstante, die wie folgt definiert ist:

(2)  G = 6{,}67408 \cdot 10^{-11}\,\mathrm{\frac{m^3}{kg \cdot s^2}}

Schon lange, bevor das newtonsche Gravitationsgesetz seinen heutigen Namen bekam, ahnte Newton höchst persönlich, dass es für solche Szenarien eigentlich nicht geeignet ist. Leider konnten seine weiterführenden Gedanken zu diesem Thema weder ihn selbst noch die Nachwelt überzeugen, so dass erst mit einigem zeitlichen Abstand Physiker wie Ernst Mach und Albert Einstein wieder Schwung in diese Problematik brachten. Das newtonsche Gravitationsgesetz bildet aus heutiger Sicht nur Grenzfälle zuverlässig ab, die nicht zeitabhängig sind und bei denen das Gravitationsfeld relativ schwach ist. Beides trifft in meiner recht einfach gehaltenen Modellphysik nicht zu. Für jedes Paar von Punktmassen 1 und 2 gilt damit zu jeder Zeit unabhängig von ihrem Abstand oder ihrer Geschwindigkeit

(3)  \vec{F}_1 = G m_1 m_2 \frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} \mathrm{\;und\;} \vec{F}_2 = G m_1 m_2 \frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3}

Damit erübrigen sich selbstverständlich Überlegungen zur Relativitätstheorie. Relativistische Masse gibt es auf diese Art nicht nicht und Punktmassen lassen sich problemlos auch weit über die Lichtgeschwindigkeit hinaus beschleunigen - ein bisschen Spaß muss sein. Auf Grund der Tatsache, dass ich nur Punktmassen simuliere, die keinen Körper haben, kann es auch nicht zu Kollisionen kommen - eine entsprechende Kollisionsabfrage entfällt demzufolge ebenfalls.

In einem Modelluniversum mit N Punktmassen gilt für Punktmasse 1 somit:

(4)  \vec{F}_1 = G m_1 \sum_{i=2}^N m_i \frac{\vec{r}_i-\vec{r}_1}{|\vec{r}_i-\vec{r}_1|^3}

Dies lässt sich in folgender Form verallgemeinern, wobei alle \vec{r} und \vec{F} zeitabhängig sind:

(5)  \vec{F}_j = G m_j \sum_{i=1}^N m_i \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3} \; \forall \; i,j \in \mathbb{Z}; 1 \leq j \leq N; i \neq j

Voranalyse für eine Simulation

Jede Punktmasse j hat drei Eigenschaften. Sie hat einen zeitabhängigen Aufenthaltsort \vec{r}_j(t) , eine zeitabhängige Geschwindigkeit \vec{v}_j(t) sowie eine konstante Masse m_j . Nach jedem diskreten Zeitschritt muss \vec{v}_j(t) aktualisiert werden, um die Punktmasse in \vec{r}_j(t) weiter zu bewegen.

Um die Geschwindigkeit neu zu berechnen muss zuerst die aus der Summe aller Kräfte auf die jeweilige Punktmasse resultierende Beschleunigung zum Zeitpunkt t berechnet werden. Zu diesem Zweck wird die letzte Gleichung durch die Ersetzung von \vec{F}_j mit m_j \vec{a_j} umgeformt:

(6)  m_j \vec{a_j} = G m_j \sum_{i=1}^N m_i \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3} \; \forall \; i,j \in \mathbb{Z}; 1 \leq j \leq N; i \neq j

Diese kann durch die Kürzung von m_j wie folgt vereinfacht werden:

(7)  \vec{a_j} = G \sum_{i=1}^N { m_i \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3} } \; \forall \; i,j \in \mathbb{Z}; 1 \leq j \leq N; i \neq j

Um eine Zeitabhängigkeit bereichert sieht diese Formel, die als erster Schritt in jeder Iteration ausgewertet werden muss, folgendermaßen aus:

(8)  \vec{a_j}(t) = G \sum_{i=1}^N { m_i \frac{\vec{r}_i(t - T)-\vec{r}_j(t - T)}{|\vec{r}_i(t - T)-\vec{r}_j(t - T)|^3} } \; \forall \; i,j \in \mathbb{Z}; 1 \leq j \leq N; i \neq j

Die obige wie auch die zwei weiteren Gleichungen bezieht sich auf den Zeitpunkt t nach einem Zeitschritt der Länge T . Im zweiten Schritt je Iteration lässt sich der neue Aufenthaltsort der Punktmasse ermitteln:

(9)  \vec{r_j}(t) = \frac{1}{2} \vec{a_j}(t) \cdot T^2 + \vec{v_j}(t - T) \cdot T + \vec{r_j}(t - T) \; \forall \; j \in \mathbb{Z}; 1 \leq j \leq N

Die neue Geschwindigkeit ergibt sich daraus im dritten und letzten Schritt je Iteration wie folgt:

(10)  \vec{v_j}(t) = \vec{a_j}(t) \cdot T + \vec{v_j}(t - T) \; \forall \; j \in \mathbb{Z}; 1 \leq j \leq N

Ein erster Entwurf

Mit diesem Wissen lässt sich an dieser Stelle schon ein Skelett für eine erste Simulation entwerfen.

 1 # Länge eines Zeitschritts
 2 T = definiere_zeitschrittlaenge()
 3 # Anzahl der zu berechnenden Zeitschritte
 4 S = definiere_zeitschrittanzahl()
 5 # Anzahl der Punktmassen
 6 N = definiere_punktmassenanzahl()
 7 # Startwerte für Punktmassen
 8 r, v, m, a = startwert(N)
 9 # Durch Zeitschritte laufen
10 for t in range(0, S):
11     # Stufe 1: Beschleunigung. Durch Punktmassen laufen
12     for j in range(0, N):
13         # Beschleunigung auf null setzen
14         a[j] = [0.0, 0.0, 0.0]
15         # Durch Punktmassen laufen
16         for i in range(0, N):
17             # Eine Masse kann sich nicht selbst anziehen
18             if j != i:
19                 # Beschleunigungen berechnen und akkumulieren
20                 a[j] = berechne_a(a[j], r[j], r[i], m[i])
21     # Stufe 2: Geschwindigkeit und Ort. Durch Punktmassen laufen
22     for j in range(0, N):
23         # Ort aktualisieren
24         r[j] = berechne_r(a[j], v[j], r[j], T)
25         # Geschwindigkeit aktualisieren
26         v[j] = berechne_v(a[j], v[j], T)
Pseudocode 1: Grober Programmentwurf in Python, der die eigentlichen Berechnungen noch bewusst außen vor lässt. Die Definition "j in range(0, N)" durchläuft alle Werte für j von 0 bis N-1.

In obigen Beispiel habe ich bewusst darauf verzichtet, eine vollständige Berechnung zu implementieren. Mein primäres Interesse gilt noch der eigentlichen Struktur der Simulation. Interessant ist, dass ich diese in zwei "Stufen" geteilt habe. Die erste Stufe übernimmt die Berechnung der Beschleunigungen. Die zweite Stufe aktualisiert in Folge dessen die Aufenthaltsorte und Geschwindigkeiten. Die Trennung der beiden Stufen ist notwendig, um nicht schon erste Punktmassen weiter zu bewegen, bevor nicht alle wirkenden Kräfte basierend auf ihren Positionen nach dem letzten Zeitschritt vollständig berechnet worden sind.

Weiterhin (negativ) auffällig ist die stark vereinfachte Implementierung von Stufe 1. Diese benötigt N (N - 1) \approx N^2 Durchläufe und entsprechend viele Aufrufe der Funktion "berechne_a". Dies führt je nach Umsetzung dieser Funktion, so fern diese ihre Ergebnisse nicht auf eine geeignete Art zwischenspeichert, zu unnötigen doppelten Berechnungen.

Eine erste strukturelle Optimierung

Bei N Punktmassen existieren N_P eindeutige Paare von Punktmassen (ohne Berücksichtigung ihrer Reihenfolge), für welche diese Beziehung bei genauer Betrachtung ausgewertet werden muss:

(11)  N_P = \frac{N (N - 1)}{2} \approx \frac{N^2}{2}

Für all diese Paare haben die Aufrufe von "berechne_a" eine Reihe an Gemeinsamkeiten, wie beispielsweise die Berechnung des Abstandes der Punktmassen. Dies lässt sich schnell zeigen. Während sich \vec{a_j} wie folgt ergibt ...

(12)  \vec{a_j} = G \cdot m_i \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}

... lässt sich die Berechnung von \vec{a_i} analog dazu mit wenigen Änderungen so ausdrücken:

(13)  \vec{a_i} = -1 \cdot G \cdot m_j \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}

Während sich die Berechnungen von \vec{a_j} und \vec{a_i} nur durch das Vorzeichen (-1) und die Masse der jeweils anderen Punktmasse unterscheiden, bleibt der rechenaufwendige Kern der Bildungsvorschriften gleich und muss je Paar nur einmal berechnet werden:

(14)  \frac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}

Eine einfache Optimierung in der Struktur des Skeletts hat somit das Potential, den Rechenaufwand von Stufe 1 von rund N^2 Durchläufen auf rund \frac{1}{2} N^2 Durchläufe und dabei auf etwas mehr als die Hälfte zu reduzieren.

 1 # Länge eines Zeitschritts
 2 T = definiere_zeitschrittlaenge()
 3 # Anzahl der zu berechnenden Zeitschritte
 4 S = definiere_zeitschrittanzahl()
 5 # Anzahl der Punktmassen
 6 N = definiere_punktmassenanzahl()
 7 # Startwerte für Punktmassen
 8 r, v, m, a = startwert(N)
 9 # Durch Zeitschritte laufen
10 for t in range(0, S):
11 	# Vorbereitung für Stufe 1
12     for j in range(0, N):
13         # Beschleunigung auf null setzen
14         a[j] = [0.0, 0.0, 0.0]
15     # Stufe 1: Beschleunigung. Durch Punktmassen laufen
16     for j in range(1, N):
17         # Durch Punktmassen laufen
18         for i in range(0, j):
19             # Beschleunigungen berechnen und akkumulieren
20             a[i], a[j] = berechne_a(a[j], r[j], m[j], a[i], r[i], m[i])
21     # Stufe 2: Geschwindigkeit und Ort. Durch Punktmassen laufen
22     for j in range(0, N):
23         # Ort aktualisieren
24         r[j] = berechne_r(a[j], v[j], r[j], T)
25         # Geschwindigkeit aktualisieren
26         v[j] = berechne_v(a[j], v[j], T)
Pseudocode 2: Grober Programmentwurf in Python mit erster struktureller Optimierung.

Indem in Stufe 1 in der äußeren Schleife j von 1 bis N und in der inneren Schleife i von 0 bis j läuft wird erreicht, dass "berechne_a" nur noch rund \frac{1}{2} N^2 mal aufgerufen wird. Ohne Verständnis für das zu berechnende physikalische Problem wäre diese Optimierung nicht denkbar gewesen. Sie an dieser Stelle durchzuführen bietet die Möglichkeit, die später zu leistenden Aufgaben von "berechne_a" richtig zu definieren - lange bevor diese Funktion mit echtem Code belebt wird.

Wie geht es weiter?

In dieser Reihe an Texten möchte ich dieses einfache Beispiel auf verschiedene Arten implementieren, dabei Stück für Stück über mehrere Größenordnungen beschleunigen und dabei entsprechende Ansätze und Techniken erklären wie auch demonstrieren. Dies wird von einfachen Schritten in unter anderem Python und Matlab bis hin zu einer vollständigen Umsetzung in C reichen, welche die Möglichkeiten moderner Hardware bis ins letzte ausreizt und sich auch auf mehrere Prozessoren und Computer verteilen lässt. Dabei ist es mir wichtig zu zeigen, wie Sie die jeweiligen Schritte gegebenenfalls auf Ihr eigenes Projekt übertragen können.

Warum ausgerechnet dieses Beispiel?

Die Formulierung des Beispiels klingt vielleicht langweilig, unrealistisch einfach oder gar fachlich hanebüchen, erzeugt aber recht eindrucksvolle Bilder, die es erlauben, das Ergebnis einer Berechnung schnell visuell zu überprüfen. Es erhebt ausdrücklich keinen Anspruch auf physikalische Korrektheit, sondern dient einzig und allein der Anschauung.

Mit einer variablen Anzahl an Punktmassen im Modelluniversum lassen sich verschiedene Aspekte der Optimierung ähnlicher Berechnungen erklären und üben. Eine geringe Anzahl an Punktmassen erlaubt es beispielsweise, einen Planeten, der um ein zentrales Gestirn kreist, zu simulieren - ein einzelnes Paar aus zwei Punktmassen, was eine schnelle prinzipielle Validierung des Codes erlaubt, ohne übermäßig viel Speicher oder Rechenleistung zu verlangen. Läuft dieser Teil einmal fehlerfrei, so lässt sich die Anzahl der simulierten Punktmassen schrittweise erhöhen, was je nach verwendeter Sprache und Implementierungsdetails immer wieder neue, grundverschiedene Engpässe hervorruft, die es mit unterschiedlichen Ansätzen zu beseitigen gilt. Das "Schöne" an diesem Beispiel ist, dass sich der Rechenaufwand mindestens quadratisch zur Anzahl der simulierten Punktmassen verhält, so dass es sich für Demonstrations- und Übungszwecke schnell und unkompliziert bis an die Grenzen handelsüblicher PCs und sogar darüber hinaus aufblasen lässt. Gleichzeitig bleibt es im Kern übersichtlich und lässt sich im Wesentlichen immer mit wenigen Zeilen Code beschreiben.

Aus mathematischer Sicht tauchen eine Reihe an typischen Operationen wenigstens einmal auf, so zum Beispiel einfache Manipulationen von Vektoren, einzelne Divisionen sowie eine Quadratwurzel. Falls gewünscht lassen sich die Berechnungen auch unter Zuhilfenahme von Winkelfunktionen durchführen, was sie zusätzlich verlangsamt und damit einen weiteren typischen Engpass einführt. Jede dieser einfachen Operationen kann, je nach Szenario und Grad der Optimierung, erhebliche Auswirkungen auf die Laufzeit der Berechnung haben. Das von mir gewählte Beispiel ist deswegen bewusst sehr einfach gehalten, um diese Effekte ohne größere Fehlersuche nach Änderungen unkompliziert studieren zu können.

Sebastian M. Ernst