Thorn Guide for the EHFinder Thorn

Original code and documentation by Peter Diener

Date

Abstract

This thorn locates the Event Horizon (EH) in an analytic or numerical spacetime 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.

1 Introduction

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 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 evolution then means that an outgoing null surface will converge exponentially to the EH. The level set function, f, is evolved according to

tf = gtiif + (gti i f)2 gtt gij i fj f gtt = βi if α2 γij i fj f, (1)

where in the second equation the lapse, shift and 3-metric has been substituted for the 4-metric. For more details on the theory and implementation see [1].

This thorn uses a level set description of the null surface, i.e.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 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 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 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 function is evolved backwards in time. The thorn can evolve more than one level set function at a time using different initial guesses for the surfaces (NOTE: this is a quite recent feature and has not yet been tested extensively). More details about the actual use of the thorn in section 7

2 Re-initialization

The evolution equation for f, equation (1), causes steepening of the gradient of f, which is undesireble numerically. For that reason, f is periodically re-initialized to a distance function. That is, the values of f are changed so that the the value of f in a grid point is equal to the (signed) distance from the grid point to the surface f = 0. This is done by evolving f according to the following evolution equation (in the parameter λ)

df dλ = f f2 + 1 |f| 1 (2)

until a steady state is achieved. This method is called the pde-method since it is basically evolving a pde. Sometimes the f = 0 surface can be moved slightly during the re-initialization procedure. This happens when 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-initialization.

There used to be another method doing the re-initialization, but it proved inferior to the pde-method and was removed. Other methods may be implemented in the future.

3 The initial shape of the surface

Currently three different choices for the initial shape of the surface are implemented. The simplest choice is a sphere in which case f is set according to

f = (x x0 )2 + (y y0 )2 + (z z0 )2 r0, (3)

where r0 is the radius of the sphere and x0, y0 and z0 are the coordinates of the center of the sphere. The second choice is a rotated and translated ellipsoid. The basic equation is here

f = x2 a2 + y2 b2 + z2 c2 1 (4)

This ellipsoid is first rotated an angle α around the z-axis, then rotated an angle β around the y-axis, then rotated an angle γ around the x-axis and finally the rotated ellipsoid is translated to have its “center” at the point (x0,y0,z0). The final possible shape of the initial surface is an ovaloid of Cassini. This was implented initially to test changing the topology in flat space. it is most likely not useful for numerical data. In this case f is set according to

f = (x2 + y2 + z2)2 + a4 2a2(x2 y2 z2) b4. (5)

There are no translation or rotations implemented for the ovaloid of Cassini.

Different initial shapes can be used for the different level set functions used in the same run.

4 Excision

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-initialized 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 ehfinder::shell_width specifies the size of the active region around f = 0. However the interior and exterior excisions are handled differently. The interior excision is most simple, since here all grid points with f < 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 < ehfinder::shell_width are in the active region. This is illustrated in Figure 1, for the case


PIC


Figure 1: Illustration of the excision regions used internally by EHFinder. The hashed regions are excised.

ehfinder::shell_width = 4, where the excised regions are hashed.

Changes to the excision regions are only done after re-initialization, since it is only at this time that f is a distance function. The excision regions can move across the grid, tracking the surface f = 0.

Also if the numerical run was done with excision using SpaceMask, the EHFinder excision region is guaranteed to completely cover the numerical excision region. The EHFinder currently only supports the old style excision mask but the support of the new style excision mask should be trivial and fast to implement.

5 Upwinding

All finite differences used in the evolution of the null surface are second order one sided differences. For that reason a ghost_size larger or equal to 2 should always be used. It is possible to choose between three different upwinding schemes. This is done by setting the parameter ehfinder::upwind_type to either intrinsic, shift or characteristic.

The intrinsic scheme, looks at the values of f itself, to determine the direction of the stencil. This is basically in order to be able to handle situations like the one illustrated in 1D in Figure 2.


PIC


Figure 2: An illustration of how to choose the upwinding stencil when f is not differntiable everywhere

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 may cause problems. The algorithm 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.

However this scheme only works in a few simple cases. For more general cases it is important to use characteristic information in order to ensure that the numerical stencil contains the domain of dependence. Therefore the characteristic scheme was introduced. Fortunately, by using the information contained in the level set functions the characteristic of the level set function can be calculated using

dxi dt = βi + α2γij jf α2 γkl k fl f. (6)

The method estimates the characteristic direction using centered finite differences in equation (6) and then recalculates the one sided finite differences in the appropriate direction.

It might happen that the upwinding direction based on the characteristic direction results in the stencil to consist of the points 1, 2, 3’ in Figure 2. However, if the re-initialization is done often enough, this turns out not to cause any problems.

The shift scheme was implemented as an improvement to the intrinsic scheme but it turned out that the characteristic scheme was superior. Therefore shift upwinding should not be used. It might be completely removed in future revisions.

For the re-initialization the default is to use the intrinsic second order scheme (the re-initialization doesn’t depend on where the surface is moving), so characteristic upwinding is not applicable. It is possible to use a first order 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.

6 The most important parameters

Here the most important parameters are described.

EHFinder also extends the following parameter from admbase in order to be able to read in initial data.

7 How to use EHFinder with numerical data

In this section I will try to describe in little more detail how EHFinder can be used to find the EH in a numerical spacetime.

7.1 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 ADMBase::metric, ADMBase::lapse and 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 EHFinder is run. If the evolution was done with a conformal factor then StaticConformal::psi has to be output as well, since it is necessary in order to reconstruct the 4-metric. However it is only necessary to output this once since it is constant during the evolution. It is not necessary to output the derivatives of the conformal factor. If excision was used in the numerical run then it is also necessary to output SpaceMask::emask1 , since EHFinder has to make sure that its internal mask covers the space mask.

At present it is necessary to output all timesteps into the same file (use IO::out_timesteps_per_file = -1, which is the default). In principle both FLEXIO and HDF5 output should be supported, but only HDF5 output has been tested. Since 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 IO::out_unchunked = "yes" to write unchunked data). 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 CactusPUGHIO/IOHDF5 for details on this). If hyperslabbing is used it is definitely necessary to do the output unchunked or recombine afterwards. An example parameter file can be seen in the par/Misner_2.2_80_3D.par.

In principle 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 EHFinder is run with the correct grid spacing and time step.

7.2 Tips for parameter choices

EHFinder is still under development and testing and can as yet 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 use at least two different initial surfaces. One surface completely inside the EH and one surface completely outside the EH. In practice it is often a good idea to use more than two different initial surfaces, since then the initial surfaces closest to the EH can be identified. Using the feature of evolving multiple level set functions at a time, this can be done while reading in the numerical data only once. 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).