BA Elias Saarmann
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.)
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 AnzahlUSER_DEF_PARAMETERSentsprechend erhöhen- den Wert von
M_centrin 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_BEGdurch(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()oderInitDomain()in den Geisterzellen (Zellen mitil<IBEGbzw.il>IEND) setzen - Wahl
userdefin pluto.ini, aber inUserDefBoundary()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,Tgasunderad, 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
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 = 1sind die Intervalle inx2-gridundx3-gridnicht irrelevant, fürGEOMETRY = SPHERICALdaher stets0 1 u 3.141592653589793bzw.0 1 u 6.283185307179586.