Nyx

Zum Beginn meiner Arbeit benutzte ich zunächst den kosmologischen Hydrodyna-mik-Code ENZO(O’Shea u. a., 2004; Bryan u. a., 2014). Aus der Vielzahl verschiede-ner kosmologischer Hydrodynamik-Codes fiel die Wahl deswegen auf ENZO, weil sein Rechengitter aus zusammenhängenden Blöcken aufgebaut ist. Jeder dieser Blöcke beinhaltet eine Vielzahl von gleichartigen, kubischen Volumenelementen³. Dies ist vor Allem im Hinblick auf die Anwendung des SGS-Modells von Schmidt u. Feder-rath (2011, fortan SF11) von Vorteil. Eine Implementierung dieses Modells wäre in baumartig strukturierten Gittern, wie in RAMSES (Teyssier, 2002), sogenannten

’moving-mesh’-Codes ohne streng reguläre Gitter wie AREPO (Springel, 2010) oder gar partikelbasierten Codes wie GADGET(Springel u. a., 2001) zumindest kaum prak-tikabel, da entspechend der Formulierung dieses Modells eine feste Auflösungsskala vorteilhaft ist, was in ENZO zumindest blockweise erfüllt ist. Leider stellte sich aber heraus, dass der in ENZO verwandte Hydrolöser nicht robust und stabil genug war, um die Problemstellung einer isolierten Scheibengalaxie zu bewältigen. So wandte ich mich der Verwendung von Nyx(Almgren u. a., 2013) zu, da auch hier blockwei-se strukturierte Rechengitter genutzt werden. Abgeblockwei-sehen von der übersichtlicheren Programmstruktur verfügt Nyxüber einen Hydrolöser, der den Anforderungen mei-ner Problemstellung besser gewachsen ist.

Nyxist ein sogenannter Hybridcode, der sowohl die strömungsmechanischen Euler-Gleichungen (siehe Kapitel 2.1.4) von baryonischem Gas auf einem numerischen Rechengitter löst, als auch andere kollisionsfreie Komponenten – im Falle kosmolo-gischer Simulationen dunkle Materie und Sterne – in Form eines N-Körpersystems von massebehafteten Partikeln berücksichtigt. Nyx basiert auf dem Softwarepaket BoxLib, das ebenfalls am Lawrence National Laboratory in Berkeley (USA) unter-halten und entwickelt wird.

Hydrodynamik

Die Gasdynamik wird in Nyx mittels eines Godunov-Verfahrens zur Zeitintegration der Feldgleichungen auf dem (adaptiven) Rechengitter berücksichtigt. Dass die hy-drodynamischen, integralen Erhaltungsgrößen unter Verwendung dieses Verfahrens tatsächlich erhalten bleiben, wird schon allein durch seine Konstruktion sicherge-stellt. Es läuft schrittweise ab, wobei für jeden Zeitschritt, dessen Länge ein Bruchteil

3Dabei kann die Größe der Volumenelemente verschiedener Gitterblöcke unterschiedlich sein, wenn das Gitter adaptiv verfeinert wird (siehe Abschnitt über AMR weiter unten).

5.1. NUMERISCHE METHODEN 49

der durch das Courant-Lewi-Strauss-Kriterium vorgegebenen kritischen Länge ist, der folgende Zyklus durchlaufen wird. Zunächst werden mittels einer konservativen Rekonstruktion zweiter Ordnung (PPM; piecewise parabolic method) der verschie-denen hydrodynamischen Felder aus den zellenzentrierten Stützdaten die Feldwerte an den jeweiligen Zellgrenzen ermittelt und dann die sich ergebenden Riemannpro-bleme an allen Grenzflächen und für alle Raumrichtungen simultan mittels eines approximativen Lösers (oder eines exakten, falls der erstere versagt) gelöst, um die Flüsse der hydrodynamischen Erhaltungsgrößen über die einzelnen Übergangsflä-chen zwisÜbergangsflä-chen den Zellen zu erhalten. Aus den Flüssen werden die entspreÜbergangsflä-chenden Inkremente für die anschließende Aktualisierung des Rechengitters berechnet. Nyx verwendet also ein numerisches Verfahren, das als ’unsplit scheme’ bezeichnet wird, also ohne getrenntes Lösen (directional splitting) der Riemannprobleme in den ein-zelnen Raumrichtungen auskommt. Die in meinen Simulationen effektiv gelösten Feldgleichungen werden im Abschnitt 5.1.3 vorgestellt.

Gravitation

Um die Eigengravitation der dynamischen Komponenten innerhalb einer Simula-tion zu berechnen, nutzt Nyx einen Poissonlöser um das Gravitationspotential zu ermitteln, woraus sich die Gravitationsbeschleunigung durch räumliche Gradienten-bildung auf dem Rechengitter bestimmen lässt. Ausgehend von der Dichte ⇢dyn des dynamischen Anteils der schweren Masse in einer Simulation wird die Poissonglei-chung

4 = 4⇡G⇢dyn (5.1)

mittels des Gauss-Seidel-Rot-Schwarz-Relaxationsverfahrens iterativ im Rahmen ei-nes Mehrgitteralgorithmus (multi grid solver) gelöst.

In meinen IDG-Simulationen wird die Eigengravitation der Gasmasse sowie der stellaren Masse auf diese Weise berücksichtigt. Die Konfiguration der dunklen Ma-terie der IDG wird als statisch angenommen. Die gravitative Wirkung des Halos dunkler Materie auf die dynamischen Komponenten wird durch ein von mir imple-mentiertes ’externes’ Beschleunigungsfeld modelliert, welches auf das interne Gravi-tationsbeschleunigungsfeld addiert wird.

Die Wirkung der Gravitationsbeschleunigung wird mittels eines Prediktor-Kor-rektor-Verfahrens auf das Gas übertragen.

Partikel

Neben Lösung der hydrodynamischen Feldgleichungen bietet Nyx die Möglichkeit auch kollisionsfreie Fluidkomponenten zu berücksichtigen, die mit sich selbst und auch mit dem Gas über die Gravitation in Wechselwirkung stehen können. Diese

werden in Form eines N-Körpersystems (N-body) behandelt. Hierbei wird das kolli-sionsfreie Fluid durch eine Anzahl von massebehafteten Partikeln repräsentiert, die sich innerhalb des Simulationsvolumens gemäß Trägheit und Gravitationsbeschleu-nigung bewegen. Zur Berechnung der Gravitation werden die Partikelmassen auf das zugrunde liegende Rechengitter mithilfe eines ’cloud in cell’ (CIC) Verfahrens projiziert, wodurch ein dem der Gasdichte vergleichbares Dichtefeld erzeugt wird.

Die aus der Lösung der Poissongleichung abgeleiteten Gravitationsbeschleunigungen werden zur Anwendung auf ein Partikel von gittergestützten Werten an die jeweilige Partikelposition interpoliert. Dieses numerische Verfahren, Partikel über Gravitati-on mittels eines Rechengitters zu koppeln, wird als PM-Methode (particle mesh method) bezeichnet. Zur Propagation der Partikel kommt in Nyx ein sogenannter

’Kick-Drift-Kick’-Algorithmus – ein spezielles Prediktor-Korrektor-Verfahren – zur Anwendung. Mittels dieses Algorithmus werden die Bewegungsgleichungen der Par-tikel schrittweise zeitlich integriert. Innerhalb eines Zeitschrittes werden zuerst die Gravitationsbeschleunigungen an den jeweiligen Positionen der Partikel ausgehend von gitterbasierten Daten interpoliert. Mithilfe der Gravitationsbeschleunigungen werden die Partikelgeschwindigkeiten zunächst um einen halben Zeitschritt inkre-mentiert (u´ = u+g t/2), um dann mittels der resultierenden Geschwindigkeiten die Partikelpositionen um einen ganzen Zeitschritt zu verschieben (´x=x+u t). Si-multan dazu wird gegebenenfalls das Gas um tpropagiert. Zuletzt, nach Lösen der Gravitation für den neuen Zustand des gesamten simulierten Systems, werden die Partikelgeschwindigkeiten mit den neuen, ebenfalls auf die Positionen interpolierten Beschleunigungen um den verbleibenden halben Zeitschritt fortentwickelt.

In kosmologischen Simulationen wird dunkle Materie in Form von Partikeln be-rücksichtigt. In meinen Simulationen isolierter Scheibengalaxien jedoch nehme ich einen statischen Halo aus dunkler Materie an, sodass nur stellare Materie als dy-namische, kollisionsfreie Komponente verbleibt, für deren Modellierung allerdings zusätzliche Programmroutinen erforderlich sind (siehe hierzu Abschnitt 5.1.3).

AMR

Da Nyx wie die numerischen Codes Maestro und Castro vom Basissoftwarepaket BoxLib abgeleitet ist, verfügtNyx über adaptive Gitterverfeinerung (AMR; adapti-ve mesh refinement). Dies ermöglicht es, Simulationen mit hoher räumlicher Auflö-sung bei gleichzeitig relativ geringem Aufwand an Rechnerressourcen durchzuführen.

Hierbei werden mittels vordefinierter Kriterien Bereiche eines Rechengitters ausge-wählt, über denen nachfolgend Gitterblöcke eingefügt werden, die über eine meist um den Faktor f_ref = 2 höhere numerische Auflösung verfügen. Diese werden mit vom gröberen Gitter interpolierten oder bereits vorhandenen feinen Daten gefüllt. Dieses Verfahren kann mehrfach angewendet werden, sodass eine Hierarchie von Gittern

5.1. NUMERISCHE METHODEN 51

mit unterschiedlichen, immer feineren Auflösungen entsteht. Im Programmablauf von Nyx werden auf diese Weise in regelmäßigen Abständen nach einer festgelegten Anzahl von Zeitschritten immer wieder aktualisierte Gitterstrukturen erzeugt, um sich der Entwicklung der hydrodynamischen Felder anzupassen. Damit die Felder in den Gittern verschiedener Auflösungen konsistent bleiben, werden die Gitterberei-che, über denen feinere Gitter liegen, mit Daten gefüllt, die sich aus der räumlichen Mittelung der darüber liegenden feineren Daten ergeben, wobei gegebenenfalls auch die Dichten der Impulse, der inneren Energie und eventuell auch der turbulenten SGS-Energie angepasst werden. Des weiteren müssen die Felder in an verfeinerte Bereiche angrenzenden Zellen auf die verfeinerten Flüsse über die Verfeinerungs-grenze hinweg korrigiert werden. Anstatt einen globalen Zeitschritt zu verwenden, der typischerweise durch die feinsten Gitterzellen vorgegeben wird, ermöglicht die-se Form von AMR, die als blockweidie-se strukturiert bezeichnet wird, sogenanntes

’sub-cycling’. Dabei können grobe Gitter mit einem relativ langen Zeitschritt fort-entwickelt werden, während feinere Gitter entsprechend ihres Verfeinerungsfaktors mehrere Zeitschrittzyklen durchlaufen, was zusätzlich den Rechenaufwand reduziert.

AMR ist besonders nützlich für Simulationen, in denen man abhängig von der Problemstellung einen großen Dynamikbereich räumlicher Skalen abdecken muss, aber die volle Auflösung nicht an allen Stellen benötigt wird, da einige Prozesse, für die möglichst hohe Auflösung notwendig ist, nur an wenigen Orten innerhalb des Simulationsvolumens relevant sind. Das spezielle Verfeinerungskriterium, das zur Anwendung kommt, muss also geeignet sein, die Bereiche zu erfassen, die für die simulierten Prozesse relevant sind. Im Falle der isolierten Scheibengalaxien habe ich ein Kriterium gewählt, das auf der Gasdichte beruht. Dieses stellt sicher, dass überall dort, wo eine Dichte von etwa 0,1 Teilchen procm³ überschritten wird, mit maximaler Auflösung gerechnet wird. Oberhalb dieser Dichte kann das ISM als das mit dem BS12-Modell beschriebene Mehrphasensystem angenommen werden, in dem auch Sternentstehung stattfinden kann.

Zusätzliche Quellterme

Abgesehen von der Behandlung rein adiabatischer Hydrodynamik und Gravitation auf dem Rechengitter bietetNyxdie Möglichkeit, weitere lokale physikalische Prozes-se wie etwa radiatives Kühlen des GaProzes-ses oder Sternentstehung zu modellieren. DieProzes-se werden als zusätzliche Quellen- und Senkenterme vor und nach dem Hydrozeitschritt auf Grundlage des Zustandes des Gases berechnet oder können als externe Quellen eingespeist werden. Diese Quellterme beeinflussen wiederum den Zustand des Gases, was mittels eines Prediktor-Korrektor-Verfahrens numerisch bewerkstelligt wird.

Auf diese Art werden in den IDG-Simulationen sowohl das Modell für nicht auf-gelöste Turbulenz nach SF11 als auch mein Modell für das turbulente, sternbildende

ISM eingebunden.