\documentclass[12pt,german]{article}
\usepackage[german]{babel}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{amsmath} % need for subequations
\usepackage{graphicx} % need for figures
\usepackage{hyperref} % use for hypertext links, including those to external documents and URLs
\usepackage{rotating}
% \setlength{\parindent}{0pt}
% \setlength{\parskip}{10pt}
\begin{document}
\title{Exercise 5\\Supernova explosion}
\author{KENN Michael, 8725258}
\date\today
\maketitle
\begin{abstract}
This report summarizes the implementation of a supernova simulation.
The simulation is based on the Euler equations, using an upwind scheme for integration.
Viscosity as well as cooling are considered.
Integration is improved by various methods such as staggered grid and van Leer integration.
Star formation is implemented as proposed in the exercise.
\end{abstract}
\section{Introduction}
This report should not be considered as an instruction for implementing a SN in hydrocode.
It is simply a summary of how I resolved any issues encountered while writing the code.
\section{Setup}
\subsection{Parametrization}
\paragraph{Normalization}:\\
Physical quantities are converted into the following units:
\begin{eqnarray*}
\left[L\right] &=& 1\,\text{pc} \\
\left[T\right] &=& 10000\,\text{years} = t_4 \\
\left[M\right] &=& M_\odot
\end{eqnarray*}
This also includes energy.
\paragraph{Physical parameters}:\\
The model is mostly based on following marginal conditions:
\begin{center}
\begin{tabular}{l|c|l}
quantity & value of quantity & value in [pc/t$_4$/M$_\odot$] \\
\hline
adiabatic index $\gamma$ & $5/3$ & 1.667 \\
supernova energy & $10^{51}\,\text{erg}$ & $5258\,M_\odot\text{pc}^2 t_4^{-2}$ \\
supernova remnant & $3\,M_\odot$ & $3\,M_\odot$ \\
inner region & $2\,\text{pc}$ & $2\,\text{pc}$ \\
molecular weight & $0.62$ & \\
ISM density & $1\,\text{cm}^{-3}$ & $0.0152\,M_\odot\text{pc}^{-3}$ \\
ISM temperate & 10000K & \\
\end{tabular}
\end{center}
\paragraph{Implementation parameters}:\\
As grid size I choose $\Delta r=0.1\,\text{pc}$.
A finer grid does not change the results significantly.
To avoid singularities, as a technical precaution, I always assume a minimal energy density $\rho\epsilon$ equivalent to 10000\,K.
This is especially crucial for the cooling.
\begin{center}
\begin{tabular}{l|c}
quantity & value \\
\hline
grid size $\Delta r$ & $0.1\,\text{pc}$ \\
minimal energy density $\rho\epsilon$ & equivalent to 10000\,K \\
viscosity factor C1 & 0.5 \\
CFL parameter for time step & 0.2 \\
\end{tabular}
\end{center}
\subsection{Implementation}
\paragraph{Variables}:\\
Besides time tracking variables, the whole status of the system at a point in time $t$ is completly described by the matrix
\begin{equation*}
\text{state}(1:4,-1:n)
\end{equation*}
The association to the physical quantities is given by the following correlations:
\begin{center}
\begin{eqnarray*}
\text{state}(1,k) & \Leftrightarrow & \rho_k^{t} \\
\text{state}(2,k) & \Leftrightarrow & m_{k+\frac{1}{2}}^{t+\frac{1}{2}} \\
\text{state}(3,k) & \Leftrightarrow & (\rho\epsilon)_{k}^{t} \\
\text{state}(4,k) & \Leftrightarrow & \rho_{k+\frac{1}{2}}^{t+\frac{1}{2}} \\
\end{eqnarray*}
\end{center}
The index $k$ represents the cell with center-distance $k\Delta r$ from the origin.
Negative indexes are required as ghost cells.
\paragraph{Physical identities}:\\
For conversion of physical quantities, the following identities are used:
\begin{eqnarray*}
T &=& (\gamma-1)\frac{\rho\epsilon}{nk_B} \\
n &=& \frac{\rho}{\mu m_P} \\
P &=& (\gamma-1)\rho\epsilon
\end{eqnarray*}
\paragraph{Gradient and Divergence}:\\
The divergence is handled as proposed in the lecture:
\begin{eqnarray*}
\nabla f(r) &=& \frac{df}{dr} \\
% \nabla\cdot\vec{a(r)} &=& \frac{1}{r^2}\frac{d(r^2 a(r))}{dr} = \frac{d(r^2 a(r))}{d(\frac{r^3}{3})}
\nabla\cdot\vec{a(r)} &=& \frac{d(r^2 a(r))}{d(\frac{r^3}{3})}
\end{eqnarray*}
\paragraph{Initialization}:\\
The overall initial SN energy is put into kinetic energy.
Nevertheless, as already mentioned, a minimum energy density equivalent to 10000\,K is maintained to ensure a realistic environment.
At discontinuities smearing over $N=5$ cells is introduced as proposed.
Since the spread of the simulation is a-priori known, allocation is done before.
\paragraph{Staggered grid}:\\
The Euler equations are forwarded consecutively as presented in the lecture on July 9th.
\begin{enumerate}
\item Euler II at ($x_{i+\frac{1}{2}}/t_{j+\frac{1}{2}}$)
\item Euler I at ($x_{i+\frac{1}{2}}/t_{j+\frac{1}{2}}$)
\item recalculate velocities at ($x_{i+\frac{1}{2}}/t_{j+\frac{1}{2}}$)
\item Euler I upwind scheme at ($x_{i}/t_{j+1}$)
\item Euler III upwind scheme at ($x_{i}/t_{j+1}$)
\end{enumerate}
\paragraph{Van Leer tangent method}:\\
The van Leer tangent method is implemented for the density equation (Euler I) and the energy equation (Euler III) as discussed in the lecture.
\paragraph{Boundary conditions}:\\
Boundary conditions are not essential for this implementation.
My approach is, to copy the most outer cell into the ghost cell.
This practice has proved to be good enough.
\paragraph{Artificial viscosity}:\\
The artifical viscosity is implemented as described in the problem description.
Pay attention, that the viscosity parameter C1 from the exercise is not dimensionless and has to be adjusted.
The time step calculation is extended accordingly.
\paragraph{Cell-shaping}:\\
The running time of this simulation is in the range of a few minutes.
Also the memory requirements are small, since we are only dealing with a one-dimensional model.
Therefore it is not required to stretch the cells in the outer regions.
However, there is far more room for performance improvements by only calculating the relevant cells around the shell.
\paragraph{Cooling}:\\
The temperate $T$ is derived from $\rho$ and $\rho\epsilon$ using the relation
\begin{eqnarray*}
T &=& \frac{\rho\epsilon}{\rho}\cdot f_T \\
f_T &\approx& 479113.1796 \\
\end{eqnarray*}
For the cooling function $\Gamma(T)$, linear interpolation between the given nodes is sufficient, since the function is smooth.
Always keep in mind
\begin{eqnarray*}
[\Gamma] &=& \frac{[M]\cdot[L]^5}{[T]^3}
\end{eqnarray*}
To get an appropriate cooling time step, I use
\begin{eqnarray*}
\left\vert\frac{\frac{d(\rho\epsilon)}{dt}\Delta t}{\rho\epsilon}\right\vert &<& 0.1 \\
\frac{d(\rho\epsilon)}{dt} &=& -n^2\Gamma(T)
\end{eqnarray*}
As an option, I also created the possibility to fall back to implicit temperature calculation by setting
\begin{eqnarray*}
\text{implicit\_cooling} &=& 1
\end{eqnarray*}
in the initial part of the code.
Implicit temperature calculation can speed up the simulation, since the cooling time step does not decrease with large temperature gradients.
My experience from this project is, that switching from explicit to implicit cooling only causes small changes in the simulation results.
\paragraph{Heating}:\\
Heating is not implemented as too much information is ambiguous to me.
\paragraph{Star formation}:\\
After 50\,Myrs the shell has traveled about 55.4\,pc.
At this point, the width of the shell is about $2\,\text{pc}$, holding a mass of about $10735\,M_\odot$.
The exercise assumes, that at a star formation rate of 10\% and a SN probability of 1\%, about $\frac{1}{1000}$ of that mass returns as energy.
The underlying energy production is $10^{51}\,\text{erg}/M_\odot$.
I inject this energy equally distributed to the shell\footnote{any other distribution as e.g. Gaussian does not show sigificant changes}.
Further I reduce the density of the shell by the star formation rate of 10\%, which is only a cosmetic correction.
The momentum within the shell has to be abandoned and therefore set to zero.
\paragraph{Output}:\\
The result of the Fortran program is input for bottom-up perl-scripts.
This is also the reason for the (on the first glimpse) not so nice format.
\section{Results}
\paragraph{Major simulations}:\\
When using artificial viscosity and cooling, after 50 Myrs, discrepancies to the theoretical predictions are less then 5\%.
It is essential to use artificial viscosity, in particular when cooling is applied.
If cooling is disregarded, the fraction between the density of the shocked gas $n_s$ and the density of the unperturbed ISM $n_0$ falls from about 5 at the end of the free expansion phase to 4.2 after 50 Myrs, coming close to the theoretical value of 4.
Cooling causes an increase in this ratio, up to a factor of 22 after about 40 Myrs, before the density of the shocked gas starts declining slowly.
The simulations with variable ISM densities illustrate, that in a denser medium the shock wave naturally expands at a slower pace.
Table \ref{TABELLE1} gives an overview of the most important results.
\begin{table}
\begin{center}
\begin{tabular}{c|c|cc|c|cc}
sim. time & n[cm$^{-3}$] & visc. & cooling & prediction & shell dist. & shell width \\
\hline
1.33\,Myrs & 1 & - & - & 14.6\,pc & 11.6\,pc & 1.4\,pc \\
1.33\,Myrs & 1 & + & - & 14.6\,pc & 14.3\,pc & 2.1\,pc \\
\hline
10\,Myrs & 1 & - & - & 36.7\,pc & 23.5\,pc & 2.1\,pc \\
10\,Myrs & 1 & + & - & 36.7\,pc & 29.1\,pc & 3.3\,pc \\
10\,Myrs & 1 & - & expl. & 36.7\,pc & 21.9\,pc & 1.3\,pc \\
10\,Myrs & 1 & + & expl. & 36.7\,pc & 26.3\,pc & 1.3\,pc \\
\hline
50\,Myrs & 1 & - & - & 58.1\,pc & 57.3\,pc & 4.9\,pc \\
50\,Myrs & 1 & + & - & 58.1\,pc & 65.9\,pc & 7.0\,pc \\
50\,Myrs & 1 & - & expl. & 58.1\,pc & 49.1\,pc & 2.0\,pc \\
50\,Myrs & 1 & + & expl. & 58.1\,pc & 55.4\,pc & 2.0\,pc \\
50\,Myrs & 1 & + & impl. & 58.1\,pc & 55.4\,pc & 2.0\,pc \\
\hline
50\,Myrs & 1 & + & - & 58.1\,pc & 65.9\,pc & 7.0\,pc \\
50\,Myrs & 2 & + & - & 48.6\,pc & 42.0\,pc & 7.7\,pc \\
50\,Myrs & 3 & + & - & 43.7\,pc & 39.3\,pc & 7.5\,pc \\
\end{tabular}
\caption{Summary of the most important simulations}\label{TABELLE1}
\end{center}
\end{table}
\paragraph{Visualization}:\\
To illustrate the extension of the shock wave for the different parameters, I created some shock wave files.
All the following seven animations can be found on my homepage in the project folder for the numerical practical\footnote{\url{http://www.kenn.at/NumericalAstronomy/}}.
\begin{itemize}
\item 10 Myrs - with viscosity - without cooling - ISM 1\,cm$^{-3}$
\item 10 Myrs - with viscosity - with cooling - ISM 1\,cm$^{-3}$
\item 50 Myrs - with viscosity - without cooling - ISM 1\,cm$^{-3}$
\item 50 Myrs - with viscosity - with and without cooling - ISM 1\,cm$^{-3}$
\item 50 Myrs - with viscosity - with cooling - ISM 1\,cm$^{-3}$ - moving grid
\item 50 Myrs - with viscosity - implicit cooling - ISM 1\,cm$^{-3}$
\item 50 Myrs - with viscosity - no cooling - ISM 1,2,3\,cm$^{-3}$
\end{itemize}
\paragraph{Star formation}:\\
I tried three different distributions to inject energy into the shell: equally distributed, triangle distributed and Gaussian distributed.
Unfortunately, in all my experiments, the system collapsed within latest 0.8 Myrs after the injection of energy due to negative densities.
I therefore also tried to spread the energy on a wider range - without sucess.
Finally I cheated by assigning stabelizing momenta.
With this approach I got physically unrealistic results.
However, the star formation process is at least theoretically covered in my source code.
\appendix
\begin{thebibliography}{5}
\bibitem{LITERATUR001}
H.~T.~Huynh,\\
{\bf Upwind Methods for the Euler Equations},\\
{\sl SIAM Journal on Numerical Analysis, Vol 32, No. 5 Oct.,1995},\\
% \url{http://}
\bibitem{LITERATUR002}
Wakana Iwakami et. al.,\\
{\bf Numerical Methods for Three-Dimensional Analysis of Shock Instability in Supernova Cores},\\
{\sl IOPscience},\\
\url{http://iopscience.iop.org/1742-6596/112/4/042021}
% \bibitem{TAG}
% Author
% {\bf Title}
% {\sl Magazin, 98-99, Januar 2006},\\
% \url{http://}
\end{thebibliography}
\end{document}