diff options
Diffstat (limited to 'doc/documentation.tex')
| -rw-r--r-- | doc/documentation.tex | 243 |
1 files changed, 204 insertions, 39 deletions
diff --git a/doc/documentation.tex b/doc/documentation.tex index 697e670..95c214a 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -81,7 +81,7 @@ \author{Original code and documentation by Peter Diener} % The title of the document (not necessarily the name of the Thorn) -\title{Thorn Guide for the {\bf EHFinder} Thorn} +\title{Thorn Guide for the {\tt EHFinder} Thorn} % the date your document was last changed, if your document is in CVS, % please us: @@ -120,16 +120,16 @@ by evolving a null surface backwards in time. The null surface is described at each time step as the 0-level iso-surface of a 3D scalar function $f$. This level set description of the surface allows, trivially, changes of the topology of the surface so it is possible to follow the merger of two (or more) black -holes into a final black hole. +holes into a final black hole. \end{abstract} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} -This thorn attempts to locate the Event Horizon (EH) in an analytic or +This thorn attempts to locate the Event Horizon (EH) in an analytic or numerical spacetime by evolving a null surface backwards in time. This -method depends on the fact that, except in cases where the coordinate are +method depends on the fact that, except in cases where the coordinates are adapted to outgoing null geodesics, an outgoing null surface started close to the EH, when evolved forward in time, is diverging exponentially from the EH. Reversing the time evolutions then means that an outgoing null surface @@ -140,22 +140,22 @@ evolved according to g^{tt}g^{ij}\partial_{i}f\partial_{j}f}}{g^{tt}}. \label{AEIDevelopment_EHFinder_evolve} \end{equation} -For more details on the theory and implementation -see~\cite{AEIDevelopment_EHFinder_Diener02}. +For more details on the theory and implementation +see~\cite{AEIDevelopment_EHFinder_Diener02}. This thorn uses a level set description of the null surface, \ie the surface -is the 0-level isosurface of a 3D scalar function, $f$, that is negative -inside and positive outside the surface. With this choice of surface +is the 0-level isosurface of a 3D scalar function, $f$, that is negative +inside and positive outside the surface. With this choice of surface description one level-set function can describe multiple surfaces at the same time, simply by having several, disconnected regions with negative values. The biggest advantage, however, is that any change of topology of -the surface is handled naturally and simply by the level-set function -changing sign. Therefore this EHFinder can follow the EH during the +the surface is handled naturally and simply by the level-set function +changing sign. Therefore this {\tt EHFinder} can follow the EH during the merger of two (or more) black holes into one black hole. To find the EH in a numerical spacetime several points have to be taken into consideration. Since the null surface has to be evolved backwards in time, the -EHFinder has to be seen as a pre-processing analysis thorn. Therefore it +{\tt EHFinder} has to be seen as a pre-processing analysis thorn. Therefore it is necessary to evolve the initial data forward in time while outputting enough 3D data, that the full 4-metric can be recovered at each timestep. The 3D data is then read back in, in reverse order, while the level-set @@ -183,13 +183,13 @@ the surface develops a narrow neck just before a topology change. For that reason, there is an option to detect when this is about to happen and undo the re-parametrization. -There is also another approximate way of performing the re-parametrization, +There is also another approximate way of performing the re-parametrization, however, even though it is faster, it is not recommended for use, since it seems to cause some instabilities and may cause larger movements of the $f=0$ surface. This method is called {\tt approx}. The two methods can be mixed, but it is not clear how well this works before -some more experiments have been done. +some more experiments have been done. \section{The initial shape of the surface} \label{AEIDevelopment_EHFinder_initial} @@ -222,9 +222,92 @@ There are no translation or rotations implemented for the ovaloid of Cassini. \section{Excision} \label{AEIDevelopment_EHFinder_excise} - +Even though the level set function, $f$, in principle can be defined +everywhere it is often advantageous to only evolve it in a certain region +around the surface $f=0$. Since $f$ is re-parametrized regularly to become +a distance function, $f$ itself can be used as a measure of the distance +from a certain grid point to the surface $f=0$. The parameter +{\tt ehfinder::shell\_width} specifies the size of the active region +around $f=0$. However the interior and exterior excision is handled +differently. The interior excision is most simply, since here all grid points +with $f<-{\mbox{\tt ehfinder::shell\_width}}$ are marked as inactive. This +should work in all cases when the excised region is everywhere concave, since +then all points on the boundary of the excised region will have enough +neighbouring active points to be able to calculate second order accurate one +sided derivatives. If the interior excised region happens to have a convex +region, this might fail. To avoid a similiar problem at the outer excised +boundary, this boundary is shaped as a rectangular box. The box is chosen +so that all points with $f<{\mbox{\tt ehfinder::shell\_width}}$ are in the +active region. This is illustrated in +Figure~\ref{AEIDevelopment_EHFinder_excisefig}, for the case +\begin{figure}[ht] + \begin{center} + \includegraphics[width=8cm]{excision.eps} + \end{center} + \caption{Illustration of the excision regions used internally by + {\tt EHFinder}. The hashed regions are excised.} + \label{AEIDevelopment_EHFinder_excisefig} +\end{figure} +{\tt ehfinder::shell\_width = 4}, where the excised regions are hashed. + +Changes to the excision regions are only done after re-parametrization, +since it is only at this time that $f$ is a distance function. The excision +regions can move across the grid, following the surface $f=0$. + +At present there is no support for using the {\tt SpaceMask} excision +mask, but this should be straightforward to implement and will be done +soon. \section{Upwinding} \label{AEIDevelopment_EHFinder_upwind} +All finite differences used in the evolution of the null surface are second +order one sided differences. For that reason a {\tt ghost\_size} larger or +equal to 2 should always be used. It is possible to choose between different +upwinding schemes depending on whether the shift is non zero or not. This is +done by setting the parameter {\tt ehfinder::upwind\_type} to either +{\tt intrinsic} (for no shift) or {\tt shift} (for non zero shift). + +The {\tt intrinsic} scheme, looks at the values of $f$ itself, to determine +the direction of the stencil. This is basically to be able to handle +situations like the one illustrated in 1D in +Figure~\ref{AEIDevelopment_EHFinder_upwindfig}. +\begin{figure}[ht] + \begin{center} + \includegraphics[width=8cm]{upwind.eps} + \end{center} + \caption{An illustration of how to choose the upwinding stencil when $f$ + is not differntiable everywhere} + \label{AEIDevelopment_EHFinder_upwindfig} +\end{figure} +If the stencil for calculating derivatives in the point labeled 1 is taken +to consist of the points 1, 2 and 3', the non differentiablility of $f$ will +cause problems. The algoithm detects this and instead uses the points 1, 2, 3 +as the stencil. This ensures that a non differentiable feature can be +maintained in the evolution. + +When there is a shift present, it can happen that the direction chosen +by the {\tt intrinsic} scheme is inconsistent with the shift direction, +causing instabilities. Since re-parametrization is done, these +instabilities are normally not allowed to develop fully, but they can +still cause distortions of the surface $f=0$ and are therefore still +damaging. For this reason it is possible to choose the upwinding direction +based solely on the shift (since $f$ is evolved backwards in time, it is +necessary to perform the upwinding in the opposite direction of the shift). + +It might happen that the upwinding direction based on the shift results +the stencil to consist of the points 1, 2, 3' in +Figure~\ref{AEIDevelopment_EHFinder_upwindfig}. This might be fixed +by increasing the frequency of re-parametrizations or by excising a +larger region inside the surface. However if these workarounds doesn't +work, it might be necessary to implement a hybrid upwinding scheme using +the shift, except at points with this problem. + +For the re-parametrization the default is to use the {\tt intrinsic} +second order scheme (the re-parametrization doesn't depend on the shift, +so {\tt shift} upwinding is not applicable. It is possible to use a +first order {\tt intrinsic} scheme, but this is, in my experience, not +accurate enough. A centered differencing scheme is also available, but +is only there for testing purposes and should never be used. These +alternative schemes will probably be removed in the future. \section{The most important parameters} Here the most important parameters are described. @@ -235,29 +318,29 @@ Here the most important parameters are described. default is {\tt normal} and should normally not be changed. This parameter may be removed in the future. \item {\tt ehfinder::eh\_metric\_type} \\ - The metric type can either be set to {\tt numerical} or {\tt analytic}. If - it is set to {\tt numerical} the EHFinder will attempt to read in the + The metric type can either be set to {\tt numerical} or {\tt analytic}. If + it is set to {\tt numerical} the {\tt EHFinder} will attempt to read in the metric from files in the directory specified by the {\tt io::recover\_dir} parameter. At present all the timesteps has to be saved in the same file. - Note that if the numerical data was produced with + Note that if the numerical data was produced with {\tt admbase::metric\_type = "static conformal"} this parameter has to be - specified again. In this case the EHFinder will attempt to also read in the - conformal factor from a file. + specified again. In this case the {\tt EHFinder} will attempt to also read + in the conformal factor from a file. If metric type is set to {\tt analytic} another thorn needs to set up the metric. It is possible to only set the metric on the initial slice, but it - is also possible to have a thorn (like thorn {\tt Exact}) set the metric at + is also possible to have a thorn (like thorn {\tt Exact}) set the metric at {\tt CCTK\_PRESTEP} if the analytic metric is time dependent. The default is at present {\tt analytic}. This should probably be changed. \item {\tt ehfinder::eh\_lapse\_type} \\ - The same for the lapse. + The same for the lapse. \item {\tt ehfinder::eh\_shift\_type} \\ - The same for the shift. + The same for the shift. \item {\tt initial\_f} \\ - The initial shape of the null surface can currently be chosen from - {\tt sphere}, {\tt ellipsoid} and {\tt cassini} as described in + The initial shape of the null surface can currently be chosen from + {\tt sphere}, {\tt ellipsoid} and {\tt cassini} as described in section~\ref{AEIDevelopment_EHFinder_initial}. The default is {\tt sphere}. \item {\tt initial\_rad} \\ - The radius of the initial sphere ($r_{0}$ in + The radius of the initial sphere ($r_{0}$ in equation~\ref{AEIDevelopment_EHFinder_sphere}). The deafault is 1. \item {\tt translate\_x} \\ How much to translate the initial surface in the $x$-direction ($x_{0}$ in @@ -287,23 +370,23 @@ Here the most important parameters are described. Rotation angle $\gamma$ for the ellipsoid around the $x$-axis. The default is 0. \item {\tt shell\_width} \\ - The width of the active evolution region. Grid points more than + The width of the active evolution region. Grid points more than {\tt shell\_width} gridspacings away from the $f=0$ surface are marked - as inactive and are not evolved as described in + as inactive and are not evolved as described in section~\ref{AEIDevelopment_EHFinder_excise}. The default is $7$ gridspacings. \item {\tt upwind\_type} \\ The type of upwinding to be used (either {\tt intrinsic} or {\tt shift}). See - the detailed description of the upwinding types in + the detailed description of the upwinding types in section~\ref{AEIDevelopment_EHFinder_upwind}. The default is {\tt intrinsic}. \item {\tt reparam\_undo} \\ - Should the re-parametrization be undone just before pinch-off or not as + Should the re-parametrization be undone just before pinch-off or not as described in section~\ref{AEIDevelopment_EHFinder_reparam}. The default is {\tt "no"}. \item {\tt re\_param\_method} \\ Choose the re-parametrization method. The choices are {\tt approx}, - {\tt pde} or {\tt mixed}. As described in - section~\ref{AEIDevelopment_EHFinder_reparam} the recommended method is + {\tt pde} or {\tt mixed}. As described in + section~\ref{AEIDevelopment_EHFinder_reparam} the recommended method is {\tt pde}. Only use {\tt approx} if you really know what you are doing. The default is {\tt pde}. \item {\tt re\_param\_int\_method} \\ @@ -317,7 +400,7 @@ Here the most important parameters are described. before giving up. Unless you are working at high resolution the default should be enough. The default is 800. \item {\tt pde\_differences} \\ - Choose the type of finite differencing used in the {\tt pde} + Choose the type of finite differencing used in the {\tt pde} re-parametrization. Don't ever use anything else than second order uwinding ({\tt upwind2}). The other choices ({\tt centered} and {\tt upwind}) are there only for testing purposes and might be removed. @@ -328,32 +411,114 @@ Here the most important parameters are described. what works best. The default (too high) is 100. \item {\tt reparametrize\_every\_approx} \\ How often to re-parametrize using the {\tt approx} re-parametrization. - Again it depends on the problem. Can be used together with + Again it depends on the problem. Can be used together with {\tt reparametrize\_every\_pde} if {\tt re\_param\_method} is set to {\tt mixed}. The default is 10. \item {\tt last\_iteration\_number} The last iteration number of the numerical simulation that produced the metric data. Active when {\tt eh\_metric\_type} is equal to {\tt numerical}. This is used in the code to figure out which - data set iteration number to read in from file + data set iteration number to read in from file ({\tt last\_iteration\_number-cctk\_iteration}). \end{itemize} -EHFinder also extends the following parameter from {\tt admbase} in order +{\tt EHFinder} also extends the following parameter from {\tt admbase} in order to be able to read in initial data. \begin{itemize} \item {\tt admbase::initial\_data} is extended with {\tt "read from file"}. \item {\tt admbase::initial\_lapse} is extended with {\tt "read from file"}. \item {\tt admbase::initial\_shift} is extended with {\tt "read from file"}. \end{itemize} -\section{How to use with numerical data} +\section{How to use {\tt EHFinder} with numerical data} \label{AEIDevelopment_EHFinder_UseThorn} - +In this section I will try to describe in little more detail how {\tt EHFinder} +can be used to find the EH in a numerical spacetime. + +\subsection{Outputting numerical data} +The first thing to make sure, is that enough data is output during the +numerical run, to be able to reconstruct the full 4-metric. The +required output therefore consists of the {\tt ADMBase::metric}, +{\tt ADMBase::lapse} and {\tt ADMBase::shift}. However if the evolution was +done with zero shift and/or lapse equal to one, it is not necessary to output +these grid functions, as long as storage is turned on and they are set to the +correct value initially when {\tt EHFinder} is run. If the evolution was done +with a conformal factor then {\tt StaticConformal::psi} has to be output as +well, since it is necessary in order to reconstruct the 4-metric. It is not +necessary to output the derivatives of the conformal factor. + +At present it is necessary to output all timesteps into the same file (use +{\tt IO::out\_timesteps\_per\_file = -1}, which is the default). In principle +both {\tt FLEXIO} and {\tt HDF5} output should be supported, but only +HDF5 output has been tested. Since {\tt EHFinder} can be run on a lot less +processors compared to the spacetime evolution, it is often advantageous to +either do unchuncked output or to recombine the output files, since it is then +possible to read the data onto a smaller number of processors (use +{\tt IO::out\_unchunked = "yes"}). If the numerical +run is larger than the EH containing region (hopefully that is the case; +otherwise the boundaries are definitely to close in), it is possible to +use hyperslabbing to just output the EH containing region (see for example +{\tt CactusPUGHIO/IOHDF5} for details on this). If hyperslabbing is used it +is definitely necessary to do the output unchunked. An example parameter +file can be seen in the {\tt par/Misner\_2.2\_80\_3D.par}. + +In principle {\tt EHFinder} should also work for downsampled (in both space +and time) data, but no experiments have been done so far to estimate the loss +of accuracy (I have always used the full resolution and done output at +every timestep). + +If hyperslabbing and/or downsampling is used, it is the users responsibility +(by specifying the right parameters in the parameter file) to ensure that +{\tt EHFinder} is run with the correct grid spacing and time step. + +\subsection{Tips for parameter choices} + +{\tt EHFinder} is still under development and testing and can as yet {\em not} +be used as a black box. But still I can give some guidelines and advice on +how to proceed. + +The first concern is to setup the initial guess for the surface. Ideally one +would like to make at least two runs with {\tt EHFinder}. One run with an +initial guess for the surface completely inside the EH and one run with an +initial guess completely outside the EH. The easiest way to get an initial +surface inside the EH, is to set up the initial guess to be completely inside +the +apparent horizon (AH). To get an initial guess that is outside of the EH +is not as easy. One way is to choose a surface, that starts to contract +everywhere when evolved according to +equation~(\ref{AEIDevelopment_EHFinder_evolve}. This of course means to do +it by trial and error. Set up some initial guess evolve it for a little while, +look at 3D output to determine if the surface is contracting everywhere and +change the initial surface if necessary. + +Then comes the question of how often to do the re-parametrization and how +much to excise. These parameters depend on the numerical data. In principle, +since the re-parametrization can move the surface, one wants to do it as +rarely as possible. On the other hand, re-parametrization is necessary +in order to keep the evolution nicely controlled (by avoiding large gradients), +so a compromise has to be found. This might require some experimentation. +Because movements of the surface during re-parametrization, usually only +occurs close to moments of topology change, it might be necessary to evolve +all the way beyond the change of topology and look at 3D output to see if +any problems occured. How often to do the re-parametrization also depends +of the width of the active region. If the active region around the surface +is narrow, it might be necessary to re-parametrize more often, since in this +case the boundaries of the active region is closer to the surface. At the +boundaries the stencil direction is dictated by the geometry and not +$f$ itself or the shift, which might cause instabilities if it is not +re-parametrized. Good initial guesses for +{\tt ehfinder::reparametrize\_every\_pde} seems to in the range 5--10, however +I have seen cases where less were required and others where more were +possible. For {\tt ehfinder::shell\_width} I normally use at least 7. + +This documentation will be updated, as input comes in from users. + +Happy event horizon finding. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{thebibliography}{9} \bibitem{AEIDevelopment_EHFinder_Diener02} - {P. Diener, {\em AEI preprint}, (2002), - 1--16. {\tt http://www.nowhere.com/}} + {P. Diener, {\em In preparation}, (2002)} +% , +% 1--16. {\tt http://www.nowhere.com/}} \end{thebibliography} % Do not delete next line |