User: rhaas Date: 2012/03/07 03:25 AM

Modified: / ET.tex /examples/collapse/ radii.py

Log: add work through most of the referee comments in examples section

File Changes:

Directory: / ============

File [modified]: ET.tex Delta lines: +95 -89 =================================================================== --- ET.tex 2012-03-07 01:15:39 UTC (rev 264) +++ ET.tex 2012-03-07 09:25:16 UTC (rev 265) @@ -1023,7 +1023,7 @@ {\tt Lorene} is often used to generate initial data involving neutron stars, where the initial data problem is sufficiently hard that finding a numerical solution requires a significant amount of computer time as well as some human -oversight in order to guide the routine to the correct solution. +vversight in order to guide the routine to the correct solution. Therefore, recommended practice is to let Lorene output data into files, and then read those into the Einstein Toolkit at the beginning of a run.

@@ -2285,11 +2285,19 @@ horizons is done by \codename{QuasiLocalMeasures}. The runs were performed 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=256\mathrm{M}$. We performed runs at 3 different +\codename{ReflectionSymmetry} and \codename{RotatingSymmetry180}), centered on +the location of the black hole. The outer +boundaries were placed at $R=256\mathrm{M}$, which places them out of causal +contact with the black whole during course of the simulation. Since the +boundary conditions employed are only approximate, this ensures that the +evolution of the black hole is not influenced by the imperfect boundary +condition. +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 +resolution on the finest (coarsest) grid. Each refined grid had twice the +resolution and half the side length of the containing grid. +The runs were performed using the tapering evolution scheme in \codename{Carpet} to avoid interpolation in time during prolongation. The initial data correspond to a rotating, stationary Kerr BH perturbed by a Brill wave\cite{Brandt:1996si} and, as such, has a non-zero @@ -2298,6 +2306,14 @@

Figure~\ref{fig:kerr_waves} shows the $\ell =2, m=0$ mode of $r\Psi_4$ extracted at $R=30\mathrm{M}$, and its numerical convergence. +In this and in the following sections, a numerical quantity $\Psi$ is +said to converge with convergence order $Q$ if for a set of numerical +stepsizes $h_1$, $h_2$, $h_3$ the difference between successive +resolutions scales as +\begin{equation} + \frac{\Psi_{h_1}-\Psi_{h_2}}{\Psi_{h_2}-\Psi_{h_3}} = \frac{h_1^Q-h_2^Q}{h_2^Q-h_3^Q},. + \label{eq:convergence-factor-definition} +\end{equation} \begin{figure} \includegraphics[width=0.9\textwidth]{waves} \includegraphics[width=0.9\textwidth]{waves_conv} @@ -2464,20 +2480,21 @@ } \end{table}

-The simulation domain spans the coordinate range -$[[x_{\rm min},x_{\rm max}],[y_{\rm min},y_{\rm max}],[z_{\rm min},z_{\rm max}]] -= [[0,120],[-120,120],[0,120]]$, where we have taken advantage of -both the equatorial symmetry (implemented by the -\codename{ReflectionSymmetry} module) and the $180\degree$ rotational -symmetry around the $z$-axis, which we apply at the $x=0$ plane using the -\codename{RotatingSymmetry180} module. -\codename{Carpet} provides a hierarchy of refined grids centered at -each puncture. Here, we used $7$ levels of refinement, -where the box edge coordinate lengths are given by -$[128,32,16,8,4,2]$ in units of the total binary mass, which is set to unity. -Note that overlapping boxes are automatically redefined by -\codename{Carpet} into one unique region before the domain -decomposition takes place. +The runs were performed +with a moving box mesh refinement provided by \codename{Carpet}, using 7 levels +of refinement on a quadrant grid (symmetries provided by +\codename{ReflectionSymmetry} and \codename{RotatingSymmetry180}). The outer +radius is located at $R = 120\mathrm{M}$, where $\mathrm{M}$ is the initial +ADM mass of the binary system. We performed runs at 5 different resolutions: +$3.125\mathrm{M}\times10^{-2} (2.0\mathrm{M})$, +$2.344\mathrm{M}\times10^{-2} (1.5\mathrm{M})$, +$1.953\mathrm{M}\times10^{-2} (1.25\mathrm{M})$, +$1.563\mathrm{M}\times10^{-2} (1.0\mathrm{M})$ and +$1.172\mathrm{M}\times10^{-2} (0.75\mathrm{M})$, where the numbers refer to +the resolution on the finest (coarsest) grid. Two set of moving boxes, each +centered on one of the black holes, are used, with the finest grid having side +length of $2\mathrm{M}$ and each coarser grid having twice the size and half +the resolution of the contained finer grid.

Figure~\ref{fig:tracks_waveform} shows the two puncture tracks throughout all phases of the binary evolution, @@ -2520,27 +2537,6 @@ solution, we ran five simulations with different resolutions, and focus our analysis on the convergence of the phase and amplitude of the Weyl scalar $\Psi_4$. -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.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.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.} -%We may define the convergence factor as the order of the Richardson extrapolation, -% -%\begin{equation} -%Q(t) = \frac{||u^{h_{\rm low}}-u^{h_{\rm med}}||} - % {||u^{h_{\rm med}}-u^{h_{\rm high}}||}, -%\end{equation} -%where $||\cdot||$ refers to an appropriate norm, and -%$h_{\rm low}$, $h_{\rm med}$ and $h_{\rm high}$ correspond to grid -%spacings of low, medium and high resolutions, respectively. -%\BCM{stops here}

Here, we consider the phase $\phi(t)$ and the amplitude $A(t)$ of the Weyl scalar $\Psi_4$ at @@ -2650,14 +2646,22 @@ 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 $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. +$M=1.4\mathrm{M}_\odot$.

+The runs are performed +with fixed mesh refinement provided by \codename{Carpet}, using 5 levels +of refinement on a quadrant grid (symmetries provided by +\codename{ReflectionSymmetry} and \codename{RotatingSymmetry180}) +The outer boundaries are placed at $R=640\mathrm{M}$, and refined boxes are +centered around the star at the origin, each doubling the resolution, with sizes +of $240\mathrm{M}$, $120\mathrm{M}$, $60\mathrm{M}$ and $30\mathrm{M}$. +We perform runs at 3 different +resolutions: the low resolution is $0.500\mathrm{M} (8.0\mathrm{M})$, medium was +$0.250\mathrm{M} (4.0\mathrm{M})$ and high was $0.125\mathrm{M} (2.0\mathrm{M})$, +where the numbers refer to the resolution on the finest (coarsest) grid. +Each refined grid had twice the resolution and half the side length of the +containing grid, with the finest grid completely covering the 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 @@ -2707,7 +2711,7 @@ 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 -directional flux-split algorithm which is of second order. Depending on +high resolution shock-capturing scheme used which is of second order. Depending on solution itself, the hydrodynamics code is only of first order in particular regions, e.g., at extrema (like the center of the star), or at the stellar surface. @@ -2738,47 +2742,52 @@ The previous examples dealt either with preexisting BHs, either single or in a binary, or with a smooth singularity free spacetime, as in the case of the TOV star. The evolution codes in the toolkit are, -however, also able to handle the dynamic formation of a singularity, +however, also able to handle the dynamic formation of a BH, that is follow a neutron star collapse into a BH. As a simple example of this process, we study the collapse of a non-rotating TOV star. We create initial data as in section~\ref{sec:tov_oscillations} using $\rho_c=3.154\times10^{-3}$ and $K_{\mathrm{ID}} = 100$, $\Gamma -= 2$, yielding a star model of gravitational mass $1.67,M_\odot$, that += 2$, yielding a star model of gravitational mass $M = 1.67,M_\odot$, that is at the onset of instability. As is common in such situations~(e.g., \cite{Baiotti:2005vi}), we trigger collapse by reducing the pressure support after the initial data have been constructed by lowering the polytropic constant $K_{\mathrm{ID}}$ from its initial value to $K = 0.98 , K_{\mathrm{ID}} = 98$. To ensure that the pressure-depleted configuration remains a solution of the Einstein constraint -equations~\eref{eqn:analysis_hamiltonian_constraint} in the presence +equations in the presence of matter, we rescale the rest mass density $\rho$ such that the total -energy density $T_{nn}$ does not change: +energy density $E$ defined in~\eref{eqn:analysis_hamiltonian_constraint} +does not change: \begin{equation} \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 being triggered by numerical noise, which would not be guaranteed to -converge to a unique solution with increasing resolution. In order to -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 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 -$8\times2^\ell,M_\odot$ for $0 \le \ell \le 4$, where the resolution -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. +converge to a unique solution with increasing resolution.

+The runs are performed 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}), +centered on the star. Refined regions are centered around the star at the +origin with the finest level having a side length of $2\mathrm{M}$ and each +coarser grid (with the exception of the coarsest level) having twice the side +length and half the resolution of the contained grid. This places the two +finest levels inside of the star; they are used to resolve the high density +region during collapse. +The outer boundaries were placed at $R=204.8\mathrm{M}$. We perform runs at 4 +different resolutions: from lowest to highest the resolution are +$0.025\mathrm{M} (3.2\mathrm{M})$, +$0.0188\mathrm{M} (2.4\mathrm{M})$, +$0.0125\mathrm{M} (1.6\mathrm{M})$ and +$4.67\times10^{-3}\mathrm{M} (0.6\mathrm{M})$, +where the numbers refer to the resolution on the finest (coarsest) grid. + 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 +sections~\ref{sec:recon} and~\ref{sec:riemann} respectively) 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 @@ -2790,11 +2799,8 @@ \caption{Coordinate radius of the surface of the collapsing star and radius 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 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 + $100$ times the atmosphere density. $R$ is the mean coordinate radius + of the apparent horizon. The lower $x$-axis displays time in code units where $M_\odot=G=c=1$, and the upper $x$-axis shows the corresponding physical time using $1,M_\odot = 4.93,\mu s$.} \end{figure} @@ -2812,30 +2818,30 @@ } \end{figure} In Figure~\ref{fig:tov_collapse_radii}, we plot the approximate -coordinate size of the star as well as the circumferential radius of -the apparent horizon that eventually forms in the simulation. The +coordinate size of the star as well as the mean coordinate radius of +the apparent horizon that eventually forms in the simulation. The apparent horizon is first found at approximately the time when the -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 different. +star's coordinate radius is twice the size of the forming apparent +horizon~\cite{Thierfelder:2010dv}. At this point the majority of the +matter is inside the apparent horizon, whose size quickly grows to its +final size while the numerical spacetime aproaches the final ``trumpet'' +configuraton.

In Figure~\ref{fig:tov_collapse_H_convergence_at0}, we display the -convergence factor obtained from -\begin{equation} - \frac{H_{h_1}-H_{h_2}}{H_{h_2}-H_{h_3}} = \frac{h_1^Q-h_2^Q}{h_2^Q-h_3^Q},, - \label{eq:convergence-factor-definition} -\end{equation} -for the Hamiltonian constraint violation at the center of the collapsing star. -Here $H_{h_i}$ is the Hamiltonian constraint -violation~\eref{eqn:analysis_hamiltonian_constraint} at the center of the -star for a run with resolution $h_i$. +convergence factor for the Hamiltonian constraint +violation~\eref{eqn:analysis_hamiltonian_constraint} at the center of +the collapsing star. Up to the time when the apparent horizon forms, the convergence order is -an expected $\approx 1.5$. At later times, the singularity forming at the +an expected $\approx 1.5$. At later times, the collapse of the lapse +function at the center of the collapsing star renders a pointwise measurement of the -convergence factor at the center impossible. +convergence factor at the center impossible. Similarly the $\Gamma$-driver +shift condition used, pushes the infalling matter beyond the +innermost grid point~\cite{Thierfelder:2010dv}. Eventually the shock +capturing scheme is unable to resolve the steep gradients and becomes +dissipative, which leads to a loss of rest mass inside of the apparent +horizon. The observed horizon mass on the other hand stays constant.

- \subsection{Cosmology}\label{sec:cosmology} The Einstein Toolkit is not only designed to evolve compact-object spacetimes, but also to solve the initial-value

Directory: /examples/collapse/ ==============================

File [modified]: radii.py Delta lines: +4 -4 =================================================================== --- examples/collapse/radii.py 2012-03-07 01:15:39 UTC (rev 264) +++ examples/collapse/radii.py 2012-03-07 09:25:16 UTC (rev 265) @@ -9,7 +9,7 @@

# load data Fx, Fy = np.genfromtxt("tov_3/star_edge.asc", comments="#", usecols=(1,2), unpack=True, invalid_raise=False) -Bx, By = np.genfromtxt("tov_3/BH_diagnostics.ah1.gp", comments="#", usecols=(1,27), unpack=True) +Bx, By = np.genfromtxt("tov_3/BH_diagnostics.ah1.gp", comments="#", usecols=(1,7), unpack=True)

# Plot basics fig = plt.figure(figsize=(6, 3)) @@ -17,8 +17,8 @@ ax = fig.add_subplot(1,1,1)

# plot -ax.plot(Fx, Fy/(2*1.63), linestyle='-', color='black', label='radius of stellar surface') -ax.plot(Bx, By/(2*1.63), linestyle='--', color='blue', label='apparent horizon radius') +ax.plot(Fx, Fy/(1.63), linestyle='-', color='black', label='radius of stellar surface') +ax.plot(Bx, By/(1.63), linestyle='--', color='blue', label='apparent horizon radius')

# plot properties #ax.set_xlim(xlim) @@ -36,7 +36,7 @@ ax2.set_xlim((ax.get_xlim()[0]/M_to_ms, ax.get_xlim()[1]/M_to_ms)) ax2.xaxis.set_major_locator(mticker.MaxNLocator(7)) ax2.xaxis.set_minor_locator(mticker.MaxNLocator(14)) -ax.set_ylabel(r'$R/R_G$') +ax.set_ylabel(r'$R [M]$') ax.yaxis.set_major_locator(mticker.MaxNLocator(8)) ax.yaxis.grid(False) set_tick_sizes(ax, 8, 4)