User: bmundim
Date: 2011/08/11 11:36 AM
Modified:
/
ET.tex
Log:
Josh's tex and bib file update
File Changes:
Directory: /
============
File [modified]: ET.tex
Delta lines: +104 -23
===================================================================
--- ET.tex 2011-08-11 12:39:27 UTC (rev 104)
+++ ET.tex 2011-08-11 16:36:33 UTC (rev 105)
@@ -186,7 +186,7 @@
available utilities such as comparison and analysis tools, typically
those specifically designed for astrophysical problems. Including them within the
Cactus computational toolkit was not felt to fit within its rapidly expanding scope. This triggered
-the creation of the Einstein Toolkit \cite{einsteintoolkit:web}. While large parts of the Einstein toolkit
+the creation of the Einstein Toolkit \cite{EinsteinToolkit:web}. While large parts of the Einstein toolkit
presently do make use of the Cactus toolkit, this is not an requirement at all,
and other contributions are welcome and have been accepted.
@@ -459,7 +459,6 @@
world-wide. To date, more than 90 peer-reviewed publication and more
than 15 student theses are based on Carpet \cite{CarpetCode:web}.
-
\subsection{Simulation Factory\pagesdone{1 Erik}}
Today's supercomputers differ significantly in
@@ -1021,7 +1020,7 @@
3-metric $\tilde{\gamma}_{ij}$, the trace $K$ of the extrinsic curvature,
the trace free extrinsic curvature $A_{ij}$ and the conformal connection
functions $\tilde{\Gamma}^i$. These are defined in terms of the
-standard ADM 4-metric $g_{ij}$, 3-metric $\gamma{ij}$, and extrinsic
+standard ADM 4-metric $g_{ij}$, 3-metric $\gamma_{ij}$, and extrinsic
curvature $K_{ij}$ by
@@ -1092,7 +1091,7 @@
S_i & := & - \frac{1}{\alpha} \left( T_{0i} - \beta^j T_{ij} \right) .
\end{eqnarray}
We have introduced the notation $\partial_0 = \partial_t -
-\beta^j\partial_j$. All quantities with a tilde $\tilde{~}$ refer to
+\beta^j\partial_j$. All quantities with a tilde involve
the conformal 3-metric $\tilde{\gamma}_{ij}$, which is used to
raise and lower indices. In particular, $\tilde{D}_i$ and
$\tilde{\Gamma}^k_{ij}$ refer to the covariant derivative and the
@@ -1271,8 +1270,9 @@
or replacing either evolution scheme.
The primary hydrodynamics evolution routine in the Einstein Toolkit is
-\codename{GRHydro}, an updated version of the public \codename{Whisky}
-code \cite{Baiotti:2004wn,Hawke:2005zw,Baiotti:2010zf,Whisky:web}. It
+\codename{GRHydro}, a code derived from the public \codename{Whisky}
+code \cite{Baiotti:2004wn,Hawke:2005zw,Baiotti:2010zf,Whisky:web}
+designed primarily by researchers at AEI and their collaborators . It
includes a high resolution shock capturing (HRSC) scheme to evolve
hydrodynamic quantities, with several different reconstruction methods
and Riemann solvers, as we discuss below. In such a scheme, we define
@@ -1296,13 +1296,15 @@
(TVD), piecewise parabolic (PPM), and essentially non-oscillatory
(ENO) methods currently implemented.
\item A Riemann problem is solved at each cell face using an
- approximate solver. Currently implemented versions include HLL
- (Harten-Lax-van Leer), Roe, and Marquina solvers.
+ approximate solver. Currently implemented versions include HLLE
+ (Harten-Lax-van Leer-Einfeldt), Roe, and Marquina solvers.
\item The conserved variables are advanced a timestep, and used to
recalculate the new values of the primitive variables.
\end{itemize}
-We discuss the GRHD formalism, the stages within a timestep, and the
-other aspects of the code below.
+We discuss the GRHD formalism, the stages within a timestep, and the
+other aspects of the code below, noting that the documentation included
+in the released version is quite extensive and covers many of these topics
+in substantially more detail.
\subsubsection{Ideal general relativistic hydrodynamics (GRHD)}
The equations of ideal GR hydrodynamics evolved by \codename{GRHydro} are
@@ -1341,12 +1343,12 @@
The {\tt GRHydro} scheme is written in a first-order hyperbolic
flux-conservative evolution system for the conserved variables
-$\hat{D}$, $\hat{S}^i$, and $\hat{\tau}$ (i.e., the only derivatives of the hydro quantities are found within the divergence expressions for flux terms, and never source terms), in terms of the primitive
+$D$, $S^i$, and $\tau$ (i.e., the only derivatives of the hydro quantities are found within the divergence expressions for flux terms, and never source terms), in terms of the primitive
variables $\rho, \epsilon, v^i$,
\begin{eqnarray}
- \hat{D} &=& \sqrt{\gamma} \rho W,\nonumber\\
- \hat{S}^i &=& \sqrt{\gamma} \rho h W^{\,2} v^i,\nonumber\\
- \hat{\tau} &=& \sqrt{\gamma} \left(\rho h W^{\,2} - P\right) - D\,,
+ D &=& \sqrt{\gamma} \rho W,\label{eq:p2c1}\\
+ S^i &=& \sqrt{\gamma} \rho h W^{\,2} v^i,\label{eq:p2c2}\\
+ \tau &=& \sqrt{\gamma} \left(\rho h W^{\,2} - P\right) - D\label{eq:p2c3}\,,
\end{eqnarray}
where $ \gamma $ is the determinant of $\gamma_{ij} $.
The evolution system then becomes
@@ -1358,10 +1360,10 @@
\end{equation}
with
\begin{eqnarray}
- \mathbf{U} & = & [\hat{D}, \hat{S}_j, \hat{\tau}], \nonumber\\
+ \mathbf{U} & = & [D, S_j, \tau], \nonumber\\
\mathbf{F}^{\,i} & = & \alpha
- \left[ \hat{D} \tilde{v}^{\,i}, \hat{S}_j \tilde{v}^{\,i} + \delta^{\,i}_j P,
- \hat{\tau} \tilde{v}^{\,i} + P v^{\,i} \right]\!, \nonumber \\
+ \left[ D \tilde{v}^{\,i}, S_j \tilde{v}^{\,i} + \delta^{\,i}_j P,
+ \tau \tilde{v}^{\,i} + P v^{\,i} \right]\!, \nonumber \\
\mathbf{S} & = & \alpha
\bigg[ 0, T^{\mu \nu} \left( \frac{\partial g_{\nu j}}{\partial x^{\,\mu}} -
\Gamma^{\,\lambda}_{\mu \nu} g_{\lambda j} \right), \nonumber\\
@@ -1371,17 +1373,96 @@
\end{eqnarray}%
%
Here, $ \tilde{v}^{\,i} = v^{\,i} - \beta^i / \alpha $ and $
-\Gamma^{\,\lambda}_{\mu \nu} $ are the 4-Christoffel symbols. The
-above equations are solved in semi-discrete fashion. The spatial
-discretization is performed by means of a high-resolution
-shock-capturing (HRSC) scheme employing a second-order accurate
-finite-volume discretization. The time integration and coupling with curvature
+\Gamma^{\,\lambda}_{\mu \nu} $ are the 4-Christoffel symbols.
+ The time integration and coupling with curvature
are carried out with the Method of
-Lines~\cite{Hyman-1976-Courant-MOL-report}.
+Lines~\cite{Hyman-1976-Courant-MOL-report}. The expressions for $\mathbf{S}$ are calculated in \codename{GRHydro} by using the definition of the extrinsic curvature to avoid any time derivatives whatsoever, as discussed in detail in the code's documentation, following a suggestion by Mark Miller based on experience with the \codename{GR3D} code.
+\subsubsection{Reconstruction techniques}
+In order to calculate fluxes at cell faces, we first must calculate function values on either side of the face. In practice, reconstructing the primitive variables yields more stable and accurate evolutions than reconstructing the conservatives. In what follows, we assume a Cartesian grid, and loop over faces along each direction in turn. We define $q^L_{i+1/2}$ to be the value of a quantity $q$ on the left side of the face between $q_i\equiv q(x_i,y,z)$ and $q_{i+1}\equiv q(x_{i+1},y,z)$, where $x_i$ is the $i$th point in the $x$-direction, and $q^R_{i+1/2}$ the right side of the same face.
+For total variation diminishing (TVD) methods, we let
+\begin{equation}
+q^L_{i+1/2} = q_i+\frac{f(q_i)\Delta x}{2};~~q^R_{i+1/2} = q_{i+1}-\frac{f(q_{i+1})\Delta x}{2}
+\end{equation}
+where $f(q_i)$ is a slope-limited gradient function, typically determined by the values of $q_{i+1}-q_i$ and $q_i-q_{i-1}$, with a variety of different forms of the slope limiter available. In practice, all try to accomplish the same task of preserving monotonicity and removing the possibility of spuriously creating local extrema. Implemented methods include minmod, superbee \cite{Roe86}, and monotonized central \cite{vanLeer77}.
+The piecewise parabolic method (PPM) is a multi-step method based around a quadratic fit to nearby points interpolated to cell faces \cite{Colella:1982ee}, for which $q^L$ and $q^R$ are generally equivalent except near shocks and local extrema. The version implemented in \codename{GRHydro} includes the steepening and flattening routines described in the original PPM papers, with a simplified dissipation procedure.
+Essentially non-oscillatory (ENO) methods use a divided differences approach to achieve third-order accuracy via polynomial interpolation \cite{Harten87,Shu99}.
+
+Both ENO and PPM yield third-order accuracy for smooth monotonic functions, whereas TVD methods typically yield second-order accurate values. Regardless of the reconstruction scheme chosen, all of these methods reduce to first order near local extrema and shocks.
+
+\subsubsection{Riemann solvers}
+
+The Riemann problem involves the solution of the equation
+\begin{equation}
+\partial_t q+\partial_i f^i(q)=0\label{eq:Riemann}
+\end{equation}
+at some point $X$ representing a discontinuity between constant states. The exact solution can be quite complicated, involving five different waves with different characteristic speeds for a hydrodynamic problem (8 for GRMHD), so \codename{GRHydro} implements three different approximate solvers to promote computational efficiency. In each case, the solution takes a self-similar form $q(\xi)$, where $\xi\equiv (x-X)/t$ represents the characteristic speed from the original shock location to the point in question in space and time.
+
+The simplest method implemented is the Harten-Lax-van Leer-Einfeldt solver \cite{Harten83,Einfeldt88} (HLL or HLLE, depending on the reference), which uses a two wave approximation to calculate the evolution along the shock front. With $\xi_-$ and $\xi_+$ the most negative and most positive wavespeeds present on either side of the interface, the solution $q(\xi)$ is assumed to take the form
+\begin{equation}
+ \label{hlle1}
+ q = \left\{ \begin{array}[c]{r c l} q^L & {\rm if} & \xi
+ < \xi_- \\ q_* & {\rm if} & \xi_- < \xi < \xi_+ \\
+ q^R & {\rm if} & \xi > \xi_+, \end{array}\right.
+\end{equation}
+\noindent with the intermediate state $q_*$ given by
+\begin{equation}
+ \label{hlle2}
+ q_* = \frac{\xi_+ q^R - \xi_- q^L - f(
+ q^R) + f(q^L)}{\xi_+ - \xi_-}.
+\end{equation}
+\noindent The numerical flux along the interface
+takes the form
+\begin{equation}
+ \label{hlleflux}
+ f(q) = \frac{\widehat{\xi}_+f(q^L) -
+ \widehat{\xi}_-f(q^R) + \widehat{\xi}_+ \widehat{\xi}_-
+ (q^R - q^L)}{\widehat{\xi}_+ - \widehat{\xi}_-},
+\end{equation}
+\noindent where
+\begin{equation}
+ \label{hlle3}
+ \widehat{\xi}_- = {\rm min}(0, \xi_-), \quad \widehat{\xi}_+ =
+ {\rm max}(0, \xi_+).
+\end{equation}
+It is these flux terms that are then used to evolve the hydrodynamic quantities.
+
+The Roe solver \cite{Roe81} involves linearizing the evolution system for the hydrodynamic evolution. Eq.~\ref{eq:Riemann}, defining the Jacobian matrix $A\equiv \frac{\partial f}{\partial q}$, and working out the eigenvalues $\lambda^i$ and left and right eigenvectors, $l_i$ and $r^j$, assumed to be orthonormalized so that $l_i\cdot r^j=\delta_i^j$. Defining the characteristic variables $w_i=l_i\cdot q$, the characteristic equation becomes
+\begin{equation}
+\partial_t w+\Lambda \partial_x w=0
+\end{equation}
+with $\Lambda$ the diagonal matrix of eigenvalues. Letting $\Delta w_i\equiv w_i^L-w_i^R=l_i\cdot (q^L-q^R)$ represent the differences in the characteristic variables across the interface, the Roe flux is calculated as
+\begin{equation}
+f(q)=\frac{1}{2}\left(f(q^L)+f(q^R)-\sum |\lambda_i| \Delta w_i r^i\right)
+\end{equation}
+where the eigenvector appearing in the summed term are evaluated for the approximate Roe average flux $q_{\rm Roe}=\frac{1}{2}(q^L+q^R)$. The Marquina flux routines use a similar approach to the Roe method, but provide a more accurate treatment for supersonic flows (i.e., those for which the characteristic wave with $\xi=0$ is within a rarefaction zone) \cite{Donat96,Aloy99b}.
+
+\subsubsection{Conservative to primitive conversion}
+
+In order to invert Eqs.~\ref{eq:p2c1} -- \ref{eq:p2c3}, solving for the primitive variables based on the values of the conservative ones, \codename{GRHydro} uses a 1-dimensional Newton-Raphson approach that solves for a consistent value of the pressure. Defining the (known) undensitized conservative variables $\hat{D}\equiv D/\sqrt{\gamma} = \rho W$, $\hat{S}^i=S^i/\sqrt{\gamma} = \rho h W^2 v^i$ and $\hat{\tau}\equiv \tau/\sqrt{\gamma} = \rho h W^2-P-\hat{D}$, as well as the auxiliary quantities $Q\equiv \rho h W^2 = \hat{\tau}+\hat{D}+P$ and $\hat{S}^2 = \gamma_{ij}\hat{S}^i\hat{S}^j = (\rho h W)^2(W^2-1)$, the former of which depends on $P$ and the latter of which is known, we find that $\sqrt{Q^2-\hat{S}^2} = \rho h W$ and thus
+\begin{eqnarray}
+\rho&=&\frac{\hat{D}\sqrt{Q^2-\hat{S}^2}}{Q}\\
+W&=&\frac{Q}{\sqrt{Q^2-\hat{S}^2}}\\
+\epsilon&=&\frac{\sqrt{Q^2-\hat{S}^2}-PW-\hat{D}}{\hat{D}}.
+\end{eqnarray}
+Given the new values of $\rho$ and $\epsilon$, one may then find the residual between the pressure and $P(\rho,\epsilon)$ and perform the Newton-Raphson step, so long as the values of $\frac{\partial P}{\partial \rho}$ and $\frac{\partial P}{\partial\epsilon}$ are known.
+
+\subsubsection{Atmospheres, boundaries, and other code details}
+
+\codename{GRHydro} uses an atmosphere, or extremely-low density floor, to avoid problems with regard to sound speeds and conservative to primitive variable conversion near the edges of matter distributions. The floor density value may be chosen in either absolute (\codename{rho\_abs\_min}) or relative (\codename{rho\_rel\_min}) terms. The atmosphere is generally assumed to have a specified polytropic EOS, regardless of the EOS describing the rest of the matter within the simulation. Whenever the numerical evolution results in a grid cell where conservative to primitive variable conversion yields negative values of either $\rho$ or $\epsilon$, the cell is reassigned to the atmosphere, with zero velocity.
+
+At present, only flat boundary conditions are supported for hydrodynamic variables, since it is generally recommended that the outer boundaries of the simulation be placed far enough away so that all cells near the edge of the computational domain represent the atmosphere.
+
+\codename{GRHydro} has the ability to advect a set of passive scalars, referred to as ``tracers'', as well as the electron fraction of a fluid, under the assumption that each tracer $X$ follows the conservation law
+\begin{equation}
+\partial_t (DX)+\partial_i(\alpha \tilde{v}^iDX)=0.
+\end{equation}
+
+\todo{Christian: Temperature Evolution?}
+
\subsection{Equations of State}
\label{sec:eoss}
%%% Written by Christian Ott
