User: knarf
Date: 2012/02/29 11:39 PM
Modified:
/
ET.tex
Log:
all my work on referee reports
File Changes:
Directory: /
============
File [modified]: ET.tex
Delta lines: +73 -58
===================================================================
--- ET.tex 2012-03-01 03:40:20 UTC (rev 252)
+++ ET.tex 2012-03-01 05:39:01 UTC (rev 253)
@@ -674,6 +674,7 @@
=1$ throughout this paper, and $M_\odot = 1$ where appropriate.
\subsubsection{ADMBase}
+\label{sec:ADMBase}
The Einstein Toolkit provides code to evolve the Einstein equations
\begin{equation}
@@ -1453,6 +1454,7 @@
\codename{GR3D} code.
\subsubsection{Reconstruction techniques}
+\label{sec:recon}
In order to calculate fluxes at cell faces, we first must calculate
function values on either side of the face. In practice, reconstructing
@@ -1491,6 +1493,7 @@
to first order near local extrema and shocks.
\subsubsection{Riemann solvers}
+\label{sec:riemann}
The Riemann problem involves the solution of the equation
\begin{equation}
@@ -1624,6 +1627,7 @@
The polytropic EOS
\begin{equation}
+\label{eq:polyEOS}
P = K\rho^\Gamma\,\,,
\end{equation}
where $K$ is the polytropic constant
@@ -1772,12 +1776,13 @@
where $n^i$ is the unit outgoing normal to the 2-surface.
The module \codename{AHFinder} provides two algorithms for locating
-\ahz{s}. The minimization algorithm~\cite{Anninos:1996ez} finds the local
-minimum of $\oint_S (\Theta - \Theta_o )^2 d^2S$ corresponding to a
+\ahz{s}. The minimization algorithm~\cite{Anninos:1996ez} finds the a
+local surface S with a minimal value for
+$\oint_S (\Theta - \Theta_o )^2 d^2S$ corresponding to a
surface of constant expansion $\Theta_o$, with $\Theta_o=0$
corresponding to the \ahz{.} For time-symmetric data, the option exists
to find instead the minimum of the surface area, which in this case
-corresponds to an \ahz{.} An alternative algorithm provided by
+corresponds to an \ahz{.} An alternative algorithm is provided by
\codename{AHFinder}, the flow algorithm~\cite{Gundlach:1997us}, on
which \codename{EHFinder} is also based.
Defining a surface as a level set $f(x^i)=r-h(\theta,\phi)=0$,
@@ -1861,7 +1866,7 @@
There are also the quasi-local measures of mass and angular momentum, from
any \ahz{s} found during the spacetime. Both
\codename{AHFinderDirect} and \codename{AHFinder} output the
-corresponding mass derived from the area of the horizon $m_H =
+corresponding irreducible mass derived from the area of the horizon $m_H =
\sqrt{A/(16\pi)}$.
The module \codename{QuasiLocalMeasures} implements
@@ -2000,11 +2005,11 @@
where $x^i$ is the puncture location and $\beta^i$ is the shift. Since the
puncture location usually does not coincide with grid points, the shift is
interpolated to the location of the puncture.
-Equation~(\eref{eq:puncturetracking}) is implemented with a simple first-order
+Equation~\eref{eq:puncturetracking} is implemented with a simple first-order
Euler scheme, accurate enough for controlling the location
of the mesh refinement hierarchy.
-Another class of objects which often needs to be tracked are neutron stars.
+Another class of objects which often needs to be tracked is neutron stars.
Here is it usually sufficient to locate the position of the maximum density
and adapt AMR resolution in these regions accordingly, coupled with the
condition that this location can only move at a specifiable maximum speed.
@@ -2048,13 +2053,14 @@
discretization procedure that specifies the number of
boundary points, their location with respect to the physical boundary,
and either the grid spacing or the number of grid points spanning the
-domain. This defines the number and
+domain. This determines the number and
location of the grid points in the discrete domain. The discrete
domain may have grid points outside of the physical domain, and may
have a non-zero boundary width. This mechanism ensures that changes in
the numerical resolution do not affect the extent of the physical
domain, i.e., that the discrete domains converge to the physical
domain in the limit of infinite resolution.
+
The Einstein Toolkit provides the \codename{CoordBase} thorn that facilitates
the definition of the simulation domain independent of the actual
evolution thorn used, allowing it to be specified at run time via parameters
@@ -2161,8 +2167,8 @@
\begin{center}
\includegraphics{bbh-boxes}
\end{center}
- \caption{Nested boxes following the individual BHs in binary
- BH merger simulation (see Section~\ref{sec:bbh-example}),
+ \caption{Nested boxes following the individual BHs in an equal-mass
+ binary BH merger simulation (see Section~\ref{sec:bbh-example}),
with the location of the individual BHs found by
\codename{PunctureTracker}. The innermost three of the nine
levels of mesh refinement used in this simulation are shown. Notice the use of
@@ -2219,8 +2225,9 @@
\label{sec:1bh-example}
As a first example, we perform simulations of a single distorted rotating BH.
We use \codename{TwoPunctures} to set up initial data for a single
-puncture of mass $M_{\mathrm{bh}}=1$ and dimensionless spin parameter
-$a = S_{\mathrm{bh}}/M_{\mathrm{bh}}^2 = 0.7$. Evolution of the data is
+puncture of mass $M_{\mathrm{bh}}=1$, a dimensionless spin parameter
+$a = S_{\mathrm{bh}}/M_{\mathrm{bh}}^2 = 0.7$ and the spin axis around the $z$ axis.
+Evolution of the data is
performed by \codename{McLachlan}, apparent horizon finding by
\codename{AHFinderDirect} and gravitational wave extraction by
\codename{WeylScal4} and \codename{Multipole}. Additional analysis of the
@@ -2228,24 +2235,24 @@
with fixed mesh refinement provided by \codename{Carpet}, using 8 levels
of refinement on a quadrant grid (symmetries provided by
\codename{ReflectionSymmetry} and \codename{RotatingSymmetry180}). The outer
-boundaries were placed at $R=256M$. We performed runs at 3 different
-resolutions: the low resolution was $0.024M (3.072M)$, medium was
-$0.016M (2.048M)$ and high was $0.012M (1.536M)$, where the numbers refer to the
-resolution on the finest (coarsest) grid. The runs where performed using the tapering
+boundaries were placed at $R=256\mathrm{M}$. We performed runs at 3 different
+resolutions: the low resolution was $0.024\mathrm{M} (3.072\mathrm{M})$, medium was
+$0.016\mathrm{M} (2.048\mathrm{M})$ and high was $0.012\mathrm{M} (1.536\mathrm{M})$, where the numbers refer to the
+resolution on the finest (coarsest) grid. The runs were performed using the tapering
evolution scheme in \codename{Carpet} to avoid interpolation in
-time during prolongation. The initial data corresponds to a rotating BH
+time during prolongation. The initial data correspond to a rotating, stationary Kerr BH
perturbed by a Brill wave and, as such, has a non-zero
gravitational wave content. We evolved the BH using 4th-order finite differencing from
-$T=0M$ until it had settled down to a stationary state at $T=120M$.
+$T=0\mathrm{M}$ until it had settled down to a stationary state at $T=120\mathrm{M}$.
Figure~\ref{fig:kerr_waves} shows the $\ell =2, m=0$ mode of $r\Psi_4$
-extracted at $R=30M$, and its numerical convergence.
+extracted at $R=30\mathrm{M}$, and its numerical convergence.
\begin{figure}
\includegraphics[width=0.9\textwidth]{waves}
\includegraphics[width=0.9\textwidth]{waves_conv}
\caption{The extracted $\ell =2, m=0$ mode of $\Psi_4$
as function of time from the high resolution run (top plot). The extraction was
- done at $R=30M$. Shown is both the real (solid black curve) and the
+ done at $R=30\mathrm{M}$. Shown is both the real (solid black curve) and the
imaginary (dashed blue curve) part of the waveform. At the bottom, we
show the
difference between the medium and low resolution runs (solid black
@@ -2263,18 +2270,19 @@
waveforms while the blue (dashed) curve shows the absolute value of the
difference between the high and medium resolution waveforms in a log-plot.
The red (dotted) curve is the same as the blue (dashed) curve, except it is
-scaled for 4th order convergence. With the resolutions used here this factor is
+scaled for 4th order convergence, demonstrating fourth-order convergence.
+With the resolutions used here this factor is
$\left (0.016^4-0.024^4\right )/\left ( 0.012^4-0.016^4\right) \approx 5.94$.
Figure~\ref{fig:kerr_waves_l4} shows similar plots for the $\ell =4, m=0$ mode
-of $r\Psi_4$, again extracted at $R=30 M$.
+of $r\Psi_4$, again extracted at $R=30 \mathrm{M}$.
\begin{figure}
\includegraphics[width=0.9\textwidth]{waves_l4}
\includegraphics[width=0.9\textwidth]{waves_l4_conv}
\caption{Real part of the extracted
$\ell =4, m=0$ mode of $\Psi_4$ as function of time (top plot) for the high
(solid black curve), medium (dashed blue curve) and low (dotted red
- curve) resolution runs. The extraction was done at $R=30M$. The bottom
+ curve) resolution runs. The extraction was done at $R=30\mathrm{M}$. The bottom
plot shows for the real part of the $\ell =4, m=0$ waveforms the
difference between the medium and low resolution runs (solid black
curve), the difference between the high and medium resolution runs
@@ -2321,14 +2329,14 @@
\end{figure}
The inset shows in more detail the differences between the different
resolutions. The irreducible mass increases by about 0.3\% during the first
-$40M$ of evolution and then remains constant (within numerical error) for the
+$40\mathrm{M}$ of evolution and then remains constant (within numerical error) for the
remainder of the evolution. The bottom plot shows the convergence of the
irreducible mass by the difference between the medium and low resolutions
(black solid curve), the difference between the high and medium resolutions
(blue dashed curved) as well as the scaled difference between the
high and medium resolutions for fourth-order (red dotted curve) and
third-order (green dash-dotted curve). The convergence is almost perfectly
-fourth-order until $T=50M$, then better than fourth-order until $T=60M$, and
+fourth-order until $T=50\mathrm{M}$, then better than fourth-order until $T=60\mathrm{M}$, and
finally between third-order and fourth-order for the remainder of the
evolution. The lack of perfect fourth-order convergence at late times may be attributed
to non-convergent errors from the puncture propagating
@@ -2367,14 +2375,14 @@
scientific interest, we have evolved a non-spinning equal-mass
BH binary system.
The initial data represent a binary system
-in a quasi-circular orbit, with an initial separation chosen
-to be $r=6M$ so we may track the later inspiral,
-plunge, merger and ring down phases of the binary evolution.
+in a quasi-circular orbit, with an initial puncture coordinate separation chosen
+to be $r=6\mathrm{M}$ so we may track the later inspiral,
+plunge, merger and ring-down phases of the binary evolution.
Table~\ref{table:BHB_ID} provides more details about
the initial binary parameters used to generate the initial data.
The \codename{TwoPunctures} module uses these initial parameters
to solve \eref{eq:twopunc_u}, the elliptic Hamiltonian constraint for
-the regular component of the conformal factor (see section~\ref{sec:twopunctures}).
+the regular component $u$ of the conformal factor (see section~\ref{sec:twopunctures}).
The spectral solution for this example was
determined by using $[n_A,n_B,n_{\phi}]=[28,28,14]$ collocation
points, and, along with the Bowen-York analytic solution for the
@@ -2417,14 +2425,14 @@
provided by the \codename{PunctureTracker} module. In the same
plot we have recorded the intersection
of the apparent horizon $2$-surface with the $z=0$ plane
-every time interval $t=10M$ during the evolution.
-A common horizon is first observed at $t=116M$. These apparent
+every time interval $t=10\mathrm{M}$ during the evolution.
+A common horizon is first observed at $t=116\mathrm{M}$. These apparent
horizons were found by the \codename{AHFinderDirect} module and their
radius and location information stored as a $2$-surface with
spherical topology by the \codename{SphericalSurface} module.
The irreducible mass and dimensionless spin of the merged BH were
calculated by the \codename{QuasiLocalMeasures} module,
-and were found to be $0.647 M$ and $-0.243 M^{-2}$, respectively.
+and were found to be $0.647 \mathrm{M}$ and $-0.243 \mathrm{M}^{-2}$, respectively.
Two modules are necessary to perform the waveform extraction.
The first one, \codename{WeylScal4}, calculates the Weyl scalar
@@ -2436,7 +2444,7 @@
mode decomposition.
Figure~\ref{fig:tracks_waveform} shows the
real and imaginary parts of the ($l=2$, $m=2$) mode for
-$\Psi_4$ extracted on a sphere centered at the origin at $R_{\rm obs} = 60M$.
+$\Psi_4$ extracted on a sphere centered at the origin at $R_{\rm obs} = 60\mathrm{M}$.
The number of grid points on the sphere was set to be
$[n_{\theta},n_{\phi}]=[120,240]$, which yields an angular
resolution of $2.6 \times 10^{-2}$ radians,
@@ -2456,11 +2464,11 @@
The mesh spacings adopted for the coarser grid in the
AMR hierarchy for these different runs were
$\{h_{\rm low},h_{\rm med},h_{\rm medh},h_{\rm high},h_{\rm higher}\}
-=\{2.0M,1.5M,1.25M,1.0M,0.75M\}$, respectively, while
+=\{2.0\mathrm{M},1.5\mathrm{M},1.25\mathrm{M},1.0\mathrm{M},0.75\mathrm{M}\}$, respectively, while
the finer grid spacings can be easily found by dividing
them by $2^k$ for the $k$th level of mesh refinement.For this example, we set
$\{h^f_{\rm low},h^f_{\rm med},h^f_{\rm medh},h^f_{\rm high},h^f_{\rm higher}\}
-=\{3.125M,2.344M,1.953M,1.563M,1.172M\}\times 10^{-2}$ for the
+=\{3.125\mathrm{M},2.344\mathrm{M},1.953\mathrm{M},1.563\mathrm{M},1.172\mathrm{M}\}\times 10^{-2}$ for the
finest grid in the different AMR hierarchies, respectively.
%JF - we assume this in the rpevious section with an equation, and can probably skip it here.
%\BCM{Convergence factor text should be moved elsewhere -- starts here.}
@@ -2477,7 +2485,7 @@
Here, we consider the phase $\phi(t)$ and
the amplitude $A(t)$ of the Weyl scalar $\Psi_4$ at
-$R_{\rm obs}=60M$. In order to take differences between
+$R_{\rm obs}=60\mathrm{M}$. In order to take differences between
the numerical values at two different grid resolutions, we use
an $8$-th order accurate Lagrange operator to interpolate the higher-accuracy finite difference solution
into the immediately coarser grid.
@@ -2515,11 +2523,11 @@
%seems to be the one corresponding to an $8$-th order accurate finite
%difference approximation.
We also indicate on both plots the time at which the gravitational
-wave frequency reaches $\omega=0.2/M$. We follow a community standard, agreed
+wave frequency reaches $\omega=0.2/\mathrm{M}$. We follow a community standard, agreed
to over the course of the NRAR\cite{NRAR:web} collaboration, that constrains
the numerical resolution so that the accumulated phase error is not
larger than $0.05$ radians at a gravitational wave frequency of
-$\omega=0.2/M$. From the plot, we assert that the phase error between the
+$\omega=0.2/\mathrm{M}$. From the plot, we assert that the phase error between the
higher and high resolutions and the one between high and medium-high
resolutions satisfies this criterion, while the phase error between
the medium-high and medium resolutions barely satisfies the criterion; and the
@@ -2535,11 +2543,11 @@
the evolution of two punctures initially located on the $x$-axis at $x=\pm 3$.
The solid blue line represents puncture 1, and the dashed red line
puncture 2. The circular dotted green lines are the intersections of the
-apparent horizons with the $z=0$ plane plotted every $10M$ during the binary
-evolution. A common horizon appears at $t=116M$. In the right panel,
+apparent horizons with the $z=0$ plane plotted every $10\mathrm{M}$ during the binary
+evolution. A common horizon appears at $t=116\mathrm{M}$. In the right panel,
we plot the real (solid blue line) and imaginary (dotted red line)
parts of the ($l=2$,$m=2$) mode of the Weyl scalar $\Psi_4$ as extracted
-at an observer radius of $R_{\rm obs}=60M$.}
+at an observer radius of $R_{\rm obs}=60\mathrm{M}$.}
\label{fig:tracks_waveform}
\end{figure}
@@ -2553,8 +2561,8 @@
dotted brown ones, the difference between high and medium-high resolutions,
while the solid blue curves represent the difference between the higher
and high resolution runs. The dotted vertical green line
-at $t=154M$ indicates the point during the evolution at which the Weyl
-scalar frequency reaches $\omega=0.2/M$. Observe that the three highest
+at $t=154\mathrm{M}$ indicates the point during the evolution at which the Weyl
+scalar frequency reaches $\omega=0.2/\mathrm{M}$. Observe that the three highest
resolutions accumulate a phase error below the standard of $0.05$ radians
required by the NRAR collaboration. }
\label{fig:amp_phs_convergence}
@@ -2573,28 +2581,28 @@
We begin our simulations with a self-gravitating fluid
sphere, described by a polytropic equation of state. This one-dimensional
-solution is obtained by the code described in section~\ref{sec:TOVSolver}, and
+initial-data set is obtained by the code described in section~\ref{sec:TOVSolver}, and
is interpolated on the three-dimensional, computational evolution grid.
This system is then evolved using the BSSN evolution system implemented in
\codename{McLachlan} and the hydrodynamics evolution system implemented in
\codename{GRHydro}.
-For the test described here, we set up a stable TOV star described by a
-polytropic equation of state $p=K\rho^\Gamma$ with $K=100$ and $\Gamma=2$,
+For the convergence test described here, we set up a stable TOV star described by a
+polytropic equation of state~\eref{eq:polyEOS} with $K=100$ and $\Gamma=2$,
and an initial central density of $\rho_c=1.28\times10^{-3}$. This model can
be taken to represent a non-rotating NS with a mass of
$M=1.4\mathrm{M}_\odot$. The computational domain is a cube of length
$640\mathrm{M}$ with a base resolution of $2\mathrm{M}$ ($4\mathrm{M}$,
$8\mathrm{M}$) in each dimension. Four additional grids refine the region
around the star centered at the origin, each doubling the resolution, with sizes
-of $120\mathrm{M}$, $60\mathrm{M}$, $30\mathrm{M}$ and $15\mathrm{M}$,
+of $240\mathrm{M}$, $120\mathrm{M}$, $60\mathrm{M}$ and $30\mathrm{M}$,
resulting in a resolution of $0.125\mathrm{M}$ ($0.25\mathrm{M}$,
$0.5\mathrm{M}$) across the entire star.
In Figure~\ref{fig:tov_rho_max} we show the evolution of the central density of
the star over an evolution time of $1300\mathrm{M}$ ($6.5\mathrm{ms}$). The
initial spike is due to the perturbation of the solution resulting from the
-interpolation onto the evolution grid. The remaining oscillations are mainly
+mapping onto the evolution grid. The remaining oscillations are mainly
due to the interaction of the star and the artificial atmosphere and are
present during the whole evolution. Given enough evolution time, the
frequencies of these oscillations can be measured with satisfactory accuracy.
@@ -2603,7 +2611,7 @@
\label{fig:tov_rho_max}
\includegraphics[width=0.9\textwidth]{rho_max}
\caption{Evolution of the central density for the TOV star. Clearly visible is
- an initial spike, produced by the interpolation of the one-dimensional equilibrium
+ an initial spike, produced by the mapping of the one-dimensional equilibrium
solution onto the three-dimensional evolution grid. The remainder of the evolution
however, the central density evolution is dominated by continuous excitations coming
from the interaction of the stellar surface with the artificial atmosphere.}
@@ -2635,7 +2643,8 @@
Within this test it is also interesting to study the convergence behavior of
the coupled curvature and matter evolution code. One of the variables often
-used for this test is the Hamiltonian constraint violation. This violation
+used for this test is the Hamiltonian constraint
+violation~\eref{eqn:analysis_hamiltonian_constraint}. This violation
vanishes for the continuum problem, but is non-zero and resolution-dependent in
discrete simulations. The expected rate of convergence of the hydrodynamics
code lies between $1$ and $2$. It cannot be higher than $2$ due to the
@@ -2645,7 +2654,7 @@
surface.
Figure~\ref{fig:tov_ham_conv} shows the order of convergence of the Hamiltonian
-constraint violation, using the three highest-resolution runs, at the stellar
+constraint violation, using the two highest-resolution runs, at the stellar
center and a coordinate radius of $r=5\mathrm{M}$ which is about half way between the
center and the surface. The observed convergence rate for most of the
simulation time lies between $1.4$ and $1.5$ at the center, and between $1.6$ and
@@ -2689,6 +2698,7 @@
\rho' + K (\rho')^2 = \rho + K_{\mathrm{ID}} \rho^2.
\label{eqn:collapse_rho_rescaled}
\end{equation}
+
Compared to the initial configuration, this rescaled star possesses a
slightly higher central density and lower pressure. This change in
$K$ accelerates the onset of collapse that would otherwise rely on
@@ -2697,7 +2707,7 @@
resolve the star as well as to push the outer boundary far enough away (so
that the star and the numerical outer boundary are not in causal
contact during the simulation) we employ a fixed mesh refinement
-scheme. The outermost box has a radius of $R_0 = 204.8\,M_\odot$ and
+scheme. The outermost box has a half-length of $R_0 = 204.8\,M_\odot$ and
a resolution of $3.2\,M_\odot$ ($2.4\,M_\odot$, $1.6\,M_\odot$,
$0.6\,M_\odot$ for higher convergence levels). Around the star,
centered about the origin, we stack $5$ extra boxes of approximate size
@@ -2705,12 +2715,16 @@
on each
level is twice that of the surrounding level. In order to resolve the large
density gradients developing during the collapse, two more levels with radii
-$4\,M_\odot$ and $2\,M_\odot$ are placed inside the star. We use the PPM
-reconstruction method and the HLLE Riemann solver to obtain second-order
+$4\,M_\odot$ and $2\,M_\odot$ are placed inside the star.
+
+We use the PPM
+reconstruction method and the HLLE Riemann solver (see
+sections~\ref{sec:recon} and~\ref{sec:riemann} resp.) to obtain second-order
convergent results in smooth regions. Due to the presence of the
density maximum at the center of the star and the non-smooth atmosphere at the
edge of the star, we expect the observed convergence rate to be somewhat lower
than second order, but higher than first order.
+
\begin{figure}
\label{fig:tov_collapse_radii}
\includegraphics[width=0.9\textwidth]{radii}
@@ -2718,7 +2732,7 @@
of the forming
apparent horizon. The stellar surface is defined as the point where $\rho$ is
$100$ times the atmosphere density. $R$ is the circumferential radius of the
- apparent horizon and $R_g = 2\,M_\star = 2\times1.63\,M_{\mathord\odot}$. An
+ apparent horizon and $R_G = 2\,M_\star = 2\times1.63\,M_{\mathord\odot}$. An
apparent horizon forms at a time roughly equal to when the mass of the star
is enclosed in its gravitational radius, forming a black hole and causally
disconnecting the evolution in the interior from the outside spacetime. The
@@ -2745,7 +2759,8 @@
star's coordinate radius approaches its Schwarzschild radius, though
one needs to keep in mind that the Schwarzschild radius is a
circumferential radius, whereas the meaning of the coordinate radius
-in our BSSN calculation is likely somewhat different.
+in our BSSN calculation is likely different.
+
In Figure~\ref{fig:tov_collapse_H_convergence_at0}, we display the
convergence factor obtained from
\begin{equation}
@@ -2785,7 +2800,7 @@
\end{equation}
which represents the familiar Kasner spacetime for a homogeneous but
anisotropically expanding universe. In the 3+1 decomposition described
-above, this reads:
+in Section~\ref{sec:ADMBase}, this reads:
\begin{widetext}
\begin{eqnarray}
\alpha(t) &=& 1 \\
@@ -2812,7 +2827,7 @@
\section{Conclusion and Future Work}
-In this article, we described the Einstein Toolkit, a collection
+In this article, we describe the Einstein Toolkit, a collection
of freely available and easy-to-use computational codes for numerical
relativity and relativistic astrophysics. The code details and example
results present in this article represent the state of the Einstein
@@ -2883,7 +2898,7 @@
%%% Several members of the CIGR team are involved
%%% in the NSF project PetaCactus, which aims to explore these possibilities.
-Besides new additions of physics modules, existing techniques require
+Besides the addition of new physics modules, existing techniques require
improvement. One example is the need for the gauge invariant
extraction of gravitational waves from simulation spacetimes as realized
by the Cauchy Characteristic Extraction (CCE) technique recently studied
