BA Elias Saarmann

• Diss auf Overleaf

Literature & Tools

• Radiating Bondi Flows I
• PLUTO manual
• pyPLUTO for evaluating the dbl-Files
• plutoplot, the successor of pyPluto
• dbl-dump, a small Python script to obtain an ASCII table from dbl-files

Bondary Conditions

in Radiating Bondi Flows I

For the outer boundry conditions \(r \rightarrow \infty\) \[ (\tilde{\rho},\mathcal{M}, s, \tilde{L}, \tilde{E}_\text{r}) \rightarrow (1, 0, s_\infty, \tilde{L}_\infty, 1) \]

With dimensionless density \( \tilde{\rho} = \frac{\rho(r)}{\rho_\infty}\), Mach number \(\mathcal{M} = \frac{v(r)}{c_s(r)}\), dimensionless entropy \( s = \frac{S}{c_\text{V}}\) and dimensionless radiation Energy density \(\tilde{E}_\text{r} = \frac{E_\text{r}}{a_\text{r}T_\infty^4}\).

in PLUTO/belt

General independent paramters:

• Density at infinity \(\rho_\infty\)
• Temperature at infinity \(T_\infty\)
• Constant opacity \(\kappa\)
• Central mass \(m\)
• molar Mass of Gas \(M\)
• Adiabatic constant \(\gamma\)

Useful dependent parameters:

• Sound velocity:

\[c_{\text{s}, \infty} = \sqrt{\frac{T_\infty \gamma R}{M}}\]

• Bondi Radius:

\[r_\text{B} = \frac{Gm}{2c_{\text{s},\infty}^2} = \frac{Gm M}{2T_\infty \gamma R}\]

Writing theoretical free parameters in terms of belt parameters:

• Optical depth:

\[\tau_\text{B} = \kappa \rho_\infty r_\text{B} = \rho_\infty^1 T_\infty^{-1} \kappa^1 m^1 M^1 \frac{G}{2 \gamma R}\]

• dimensionless cooling time:

\[\beta = \frac{1}{4 \gamma (\gamma -1)} \frac{\rho_\infty c_{\text{s},\infty}^3}{ac T_\infty^4} \frac{1}{\tau_\text{B}} = \rho_\infty^0 T_\infty^{-\frac{3}{2}} \kappa^{-1} m^{-1} M^{-\frac{5}{2}} \frac{\gamma^{\frac{3}{2}} R^\frac{5}{2}}{4 (\gamma -1) G a c}\]

Parameters for belt:

• Central Mass:

\[m = m_\odot\]

• Optical density:

\[\kappa = \]

• Adiabatic constant/degrees of freedom

\[\gamma = \frac{7}{5} \rightarrow \text{degrees of freedom} = 5\]

• inner radius:

\[r_\text{s} = 10^{-2} r_B\]

• outer radius:

\[r_\text{s} = 10^{-2} r_B\]

• molecular mass:

\[m_\text{molec} = \]

Boundary conditions for fields in belt...

...at outer radius \(r_\text{o}\):

\[ p(r_\text{o}) = \frac{\rho_\infty k_\text{b} T_\infty}{m_\text{molec}}\]

\[ \rho (r_\text{o}) = \rho_\infty\]

\[\partial_r v(r_\text{o}) = 0\]

\[E_\text{r}(r_\text{o}) = a_\text{r} T_\infty^4\]

...at inner radius \(r_\text{s}\):

\[\partial_r p(r_\text{s}) = 0\]

\[\partial_r \rho (r_\text{s}) = 0\]

\[\partial_r v(r_\text{s}) = 0\]

\[\partial_r E_\text{r}(r_\text{s}) = 0\]

(While the zero gradient conditions do not match theoretical curves they should still be appropriate as the will only lead to local devastationdeviations near the respective boundary.)

Final Analytical Approximation

Analytical Approximation for \(f_\text{acc}\)

For \( \frac{\tilde{L}_\infty}{10 \beta} < \tau_\text{B} < \tilde{L}_\infty^{-1}\) \[ f_\text{acc} = \text{No Solution}\]

For \(10 \tau_\text{B} \beta < \tilde{L}_\infty < \text{min}(10, \tau_\text{B}^{-1})\) \[ f_\text{acc} = 1\]

For \(\tilde{L}_\infty > \text{max}(10, 100\tau_\text{B})\) \[ f_\text{acc} = \left(\frac{\tilde{L}_\infty}{10}\right)^{-\frac{5}{4}}\]

For \(\tau_\text{B}^{-1} < \tilde{L}_\infty < 100\tau_\text{B} \) and \(\tilde{L}_\infty > \tau_\text{B}^\frac{5}{11} (10\beta)^\frac{8}{11} \) \[ f_\text{acc} = (\tilde{L}_\infty \tau_\text{B})^{-\frac{5}{8}}\]

For \(\tau_\text{B}^{-1} < \tilde{L}_\infty < 100\tau_\text{B} \) and \(\tilde{L}_\infty < \tau_\text{B}^\frac{5}{11} (10\beta)^\frac{8}{11}\) \[ f_\text{acc} = (10\beta \tau_\text{B}^2)^{-\frac{5}{11}}\]

Analytical Approximation for \(\tilde{L}_\infty\)

\[\tilde{L}_\infty = \frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta \tau_\text{B} f_\text{acc} \]

For \( \frac{\tilde{L}_\infty}{10 \beta} < \tau_\text{B} < \tilde{L}_\infty^{-1}\) \[ \tilde{L}_\infty = \text{No Solution}\]

For \(10 \tau_\text{B} \beta < \tilde{L}_\infty < \text{min}(10, \tau_\text{B}^{-1})\) \[ \tilde{L}_\infty = \frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta \tau_\text{B}\]

For \(\tilde{L}_\infty > \text{max}(10, 100\tau_\text{B})\) \[ \tilde{L}_\infty = \frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta \tau_\text{B} \left(\frac{\tilde{L}_\infty}{10}\right)^{-\frac{5}{4}}\]

\[ \tilde{L}_\infty = \left(\frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta \tau_\text{B} 10^{\frac{5}{4}}\right)^\frac{4}{9}\]

For \(\tau_\text{B}^{-1} < \tilde{L}_\infty < 100\tau_\text{B} \) and \(\tilde{L}_\infty > \tau_\text{B}^\frac{5}{11} (10\beta)^\frac{8}{11} \) \[ \tilde{L}_\infty = \left(\frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta \tau_\text{B}^\frac{3}{8} \right)^{\frac{8}{13}}\]

For \(\tau_\text{B}^{-1} < \tilde{L}_\infty < 100\tau_\text{B} \) and \(\tilde{L}_\infty < \tau_\text{B}^\frac{5}{11} (10\beta)^\frac{8}{11}\) \[ \tilde{L}_\infty = \frac{56}{6}\frac{r_\text{B}}{r_\text{s}}\beta^{\frac{6}{11}} \tau_\text{B}^{\frac{1}{11}} 10^{-\frac{5}{11}}\]

das Reference-System

(Lothar)

Die Definitionen von UNIT_DENSITY, UNIT_LENGTH und UNIT_VELOCITY in user_defined_parameters.h legen die Umrechungsfaktoren von g/cm³, cm und cm/s in ein vom User gewähltes Einheitensystem (den "Code-Einheiten") fest. Daraus lassen sich dann die Umrechungsfaktoren auch für alle anderen Größen definieren, was durch den Aufruf InitializeReferenceSystem(); in init.c durchgeführt wird. Der Umrechnungsfaktor heißt dann immer Reference.... Daher ist es eine Standard-Aktion, einen User-Parameter aus pluto.ini zu nehmen (g_inputParam[...]) und ihn meist sofort durch das entsprechende Reference... zu dividieren, um ihn in Code-Einheiten zu haben.

Die ReferenceTemperature fällt ein bisschen raus, da sie auf 1 fixiert ist, unabhängig von obigen drei Basisgrößen. Die Code-Einheit für die Temperatur ist also stets Kelvin und die Division durch die ReferenceTemperature wird nur "der Konsistenz halber" gemacht.

Code-Anpassungen...

(Lothar)

für konstante Zentralmasse

• M_centr (o.Ä.) als weiteren Parameter in user_defined_parameters.h definieren und die Anzahl USER_DEF_PARAMETERS entsprechend erhöhen
• den Wert von M_centr in pluto.ini unter [Parameters] setzen (in Gramm)
• body-force.c aus belt/src/Misc in den Run-Folder kopieren und in Zeile 64 die M_X1_BEG durch (g_inputParam[M_centr]/ReferenceMass) ersetzen

für Dirichlet-RBn

Abgesehen von reflective zur Nullsetzung der Normalkomponenten von \(\vec v\) und \(\vec B\) bietet PLUTO keine Dirichlet-RBn. Workaround:

• gewünschte Randwerte in Init() oder InitDomain() in den Geisterzellen (Zellen mit il<IBEG bzw. il>IEND) setzen
• Wahl userdef in pluto.ini, aber in UserDefBoundary() am entsprechenden Rand gar keine Manipulationen vornehmen

aktuelle Probleme

Simulationsprobleme

Versuch \(\tau_\text{B} = 10^{-3}\) und \(\beta = 10^{-3}\) zu reproduzieren

• Problem: Code bricht ab, dt wird zu klein
  • Simulationsparameter: \[\rho_\infty = 3 \cdot 10^{-20} g/cm^3\] \[T_\infty = 2.7 K\] \[\kappa = 0.1 cm^2/g\] \[m = 1 m_\odot\] \[M = 2.01588 g/Mol\] \[r_\text{s} = 3 au\] \[r_\text{o} = 206265.0 au\] grid cells = 114
  • theoretische Parameter:

\[\tau_\text{B} = 1.28 \cdot 10^{-3}\] \[\beta = 1.69 \cdot 10^{-3}\] \[r_\text{B} = 28450.69 au\]

Versuch \(\tau_\text{B} = 10^{3}\) und \(\beta = 10^{-3}\) zu reproduzieren

• Problem: Starke Schwankungen nahe Innenrand für prs , Tgas und erad , Zero-Gradient Randbedingungen werden nicht eingehalten. Negative Drücke in Log File:
  • Simulationsparameter: \[\rho_\infty = 3 \cdot 10^{-14} g/cm^3\] \[T_\infty = 2.7 K\] \[\kappa = 0.1 cm^2/g\] \[m = 1 m_\odot\] \[M = 2.01588 g/Mol\] \[r_\text{s} = 3 au\] \[r_\text{o} = 206265.0 au\] grid cells = 114
  • theoretische Parameter:

\[\tau_\text{B} = 1.28 \cdot 10^{3}\] \[\beta = 1.69 \cdot 10^{-3}\] \[r_\text{B} = 28450.69 au\]

Lösungsansätze:

• grid cells auf 150 erhöhen

Simulationsschritte werden schnell sehr klein dt = 1e-9 = 10^-9 und Simulation kommt kaum voran

• Parameter varrieren, so dass Bondi Radius kleiner wird
  • Simulationsparameter: \[\rho_\infty = 3 \cdot 10^{-12} g/cm^3\] \[T_\infty = 10 K\] \[\kappa = 0.01 cm^2/g\] \[m = 1 m_\odot\] \[M = 2.01588 g/Mol\] \[r_\text{s} = 3 au\] \[r_\text{o} = 206265.0 au\] grid cells = 150
  • theoretische Parameter:

\[\tau_\text{B} = 3.45 \cdot 10^{3}\] \[\beta = 2.38 \cdot 10^{-3}\] \[r_\text{B} = 7681.7 au\]

Läuft, allerings wieder viele negative Drücke in der Log file und Ergebnisse mit Oszilation und Verletzung der Zero Gradient Bedingung bei Druck, Temperatur und Strahlungsenergiedichte.

• mehr Grid Zellen, optische Dichte \(\tau = \kappa \rho \Delta x\) der einzelnen Zellen sollte klein sein
  • Parameter wie oben, 300 Zellen -> Oszilationen werden schwächer, power law erkennbar, immernoch Probleme mit Zero-Gradient-Randbedingung
  • Parameter wie oben, 1000 Zellen -> nur noch sehr sehr schwache Oszillationen, power law klar erkennbar, Zero-Gradient scheinbar efüllt

Allerdings: Massenflux steigt linear, kein Steady State erreicht

-> Laufzeit erhöhen:

• \(t_\text{final} = 4\cdot10^6 \text{a}\) mit 300 Zellen:
  • Simulation läuft etwa 12.5 h
  • große Schwankungen in allen Daten
• \(t_\text{final} = 1.3\cdot10^6 \text{a}\) mit 1000 Zellen:
  • Simulation läuft etwa 29 h
  • keine Schwankungen
• in beiden Fällen
  • Akkretionsrate nähert sich Asymptote an bei \(f_\text{acc} \approx 1\) (entspricht \(\tilde{L}_\infty \approx 2\cdot10^5\)) während nach Theorie \(f_\text{acc} \approx 4\cdot10^{-4}\) (entspricht \(\tilde{L}_\infty \approx 8\cdot10^1\)) erwartet wird
  • in vx_1 ist Schock erkennbar der noch nicht ganz durchgelaufen ist
  • Dichteverteilung fehlt das Minimum aus Fig 6 im Paper
  • T und E_rad ähneln Verlauf in Figure 6

-> Fragen:

• Ist es Korrekt Innenrand X_Beg als \(r_\text(s)\) für Akkretionsluminosität zu behandeln?
• Kann ich die Luminosität direkt in den Daten sehen oder muss ich diese über Akkretionsrate berechnen?
  • Könnte es sein dass Probleme beim Strahlungstransport vorliegen?
  • Warum brauchen die Simulationen so lange?
  • Spielt der Schock der noch unterwegs ist eine größere Rolle oder ändert sich an \(f_\text{acc}\) später auch nichts mehr?

Versuch \(\tau_\text{B} = 10^{-1}\) und \(\beta = 10^{-3}\) zu untersuchen

• Simulationsparameter:

\[\rho_\infty = 3 \cdot 10^{-18} g/cm^3\] \[T_\infty = 2.7 K\] \[\kappa = 0.1 cm^2/g\] \[m = 1 m_\odot\] \[M = 2.01588 g/Mol\] \[r_\text{s} = 3 au\] \[r_\text{o} = 206265.0 au\] grid cells = 114

• theoretische Parameter:

\[\tau_\text{B} = 1.28 \cdot 10^{-1}\] \[\beta = 1.69 \cdot 10^{-3}\] \[r_\text{B} = 28450.69 au\]

• Problem: Zeitschritte werden sehr klein (dt = 1e-9) und die Simulation braucht ewig: nach fast 3e6 Simulationsschritten immer noch bei t = 0.23.
  • Simulationsparameter: \[\rho_\infty = 3 \cdot 10^{-18} g/cm^3\] \[T_\infty = 2.7 K\] \[\kappa = 0.1 cm^2/g\] \[m = 1 m_\odot\] \[M = 2.01588 g/Mol\] \[r_\text{s} = 3 au\] \[r_\text{o} = 206265.0 au\] grid cells = 114
  • theoretische Parameter:

\[\tau_\text{B} = 1.28 \cdot 10^{-1}\] \[\beta = 1.69 \cdot 10^{-3}\] \[r_\text{B} = 28450.69 au\]

auf Elias' Laptop (Ubuntu 26)

mpirun nötig

./belt ohne mpirun bleibt direkt "stecken".

-O2 nötig

./belt ohne -O2 compiliert bleibt nach dem Profiling-Output "stecken" oder bricht mit Fehlermeldung ab.

pluto.ini-Syntax

(gelöst)

• Alle Parameter müssen angegeben werden, auch wenn gar kein Zugriff mittels g_inputParam[] erfolgt.
• Die Reihenfolge der Parameter muss der Nummerierung in user_defined_parameters.h folgen.
• Auch für DIMENSIONS = 1 sind die Intervalle in x2-grid und x3-grid nicht irrelevant, für GEOMETRY = SPHERICAL daher stets 0 1 u 3.141592653589793 bzw. 0 1 u 6.283185307179586.