AHFinderDirect – A Fast Apparent Horizon Finder
Jonathan Thornburg <jthorn@aei.mpg.de>
$$Date$$
Thorn AHFinderDirect locates apparent horizons (or more generally, closed 2surfaces with
${S}^{2}$
topology having any desired constant expansion) in a numerically computed slice using a direct
method, posing the apparent horizon equation as an elliptic PDE on angularcoordinate space. This
is very fast and accurate, but requires a “reasonable” initial guess. This thorn guide describes how
to use the thorn.
1 Introduction
A “marginally trapped surface” is a closed 2surface in a slice whose congruence of futurepointing
outgoing null geodesics has zero expansion. There may be several such surfaces, some nested inside
others; an “apparent horizon” is an outermost marginally trapped surface. In terms of the usual
$3+1$
variables, an apparent horizon satisﬁes the equation
$$\Theta \equiv {\nabla}_{i}{n}^{i}+{K}_{ij}{n}^{i}{n}^{j}K=0$$  (1) 
where ${n}^{i}$
is the outwardpointing unit normal to the apparent horizon, and
${\nabla}_{i}$
is the covariant derivative operator associated with the 3metric
${g}_{ij}$ in the slice. (See [?] for a derivation
of equation $\left(1\right)$.) (Optionally, you
can replace the right hand side of $\left(1\right)$
by any speciﬁed nonzero constant, i.e. you can ﬁnd a surface of constant (in general nonzero) expansion.; this is
dicsussed in section 4.8.)
Thorn AHFinderDirect ﬁnds an apparent horizon by numerically solving
equation $\left(1\right)$. It requires as input
the usual Cactus 3metric ${g}_{ij}$
and extrinsic curvature ${K}_{ij}$, (and
optionally the conformal factor $\psi $
if the StaticConformal metric semantics are used), and produces as output the Cactus
$\left(x,y,z\right)$
coordinates of a large number of points on the apparent horizon, together with some auxiliary information like
the apparent horizon area and centroid position, and the irreducable mass associated with the
area.
Besides this thorn guide, the other main sources of information on AHFinderDirect are the comments in the
param.ccl ﬁle, the paper [?], and to a lesser extent the paper [?]. As a courtesy, I ask that both
these papers be cited in any published research which uses this thorn, or which uses code from this
thorn.
2 What AHFinderDirect Needs
There are some restrictions on the spacetime, or more precisely on each slice where you want to ﬁnd apparent
horizons, which are necessary in order for AHFinderDirect to work:
 AHFinderDirect requires that the Cactus geometry (${g}_{ij}$,
${K}_{ij}$,
and optionally $\psi $)
be nonsingular in a neighborhood of the apparent horizon. In particular, this means that it quite
certainly will not work for spacetimes/slicings which have a singular geometry on the horizon, such
as Schwarzschild/Schwarzschild and Kerr/BoyerLindquist.
 Less obviously, this also means that if there is a singularity in the geometry somewhere near the
apparent horizon, then you need to have a high enough Cactus 3D grid resolution that the geometry
interpolation doesn’t “see” the singularity. (If AHFinderDirect “sees” the singularity, it may
“just” fail to ﬁnd the horizon, and/or it may report that the interpolated ${g}_{ij}$
fails to be positive deﬁnite or even contains NaNs.)
 At the moment AHFinderDirect and the Cactus interpolators don’t know how to avoid an
excised region, so if the apparent horizon (or any trial horizon surface as the algorithm is iterating
towards the apparent horizon) gets too close to an excised region, you’ll get garbage results as the
interpolator tries to interpolate data from the excised region. I plan to ﬁx this sometime soon.
 AHFinderDirect requires that any apparent horizon it’s going to (try to) ﬁnd must be a
“Strahlkörper” (literally “ray body”, or more commonly “starshaped region”) relative to some local
coordinate origin (which you must specify). A Strahlkörper is deﬁned by Minkowski ([?, p. 108])
as
a region in $n$dimensional
Euclidean space containing the origin and whose surface, as seen from the origin, exhibits
only one point in any direction.
In other words, using polar spherical coordinates relative to the local
coordinate origin, the apparent horizon’s shape must be parameterizable as
$r=h\left(angle\right)$ for some singlevalued
function $h:{S}^{2}\to {\Re}^{+}$.
(AHFinderDirect uses precisely this parameterization.)
There are also some restrictions on your Cactus conﬁguration and run; here’s what works and what doesn’t:
 I strongly recommend using a currentCVS checkout of the Cactus ﬂesh and of all relevant thorns.
I haven’t tested AHFinderDirect at all with older versions of the ﬂesh or other thorns.
 AHFinderDirect works ﬁne with the PUGH unigrid driver and with the Carpet
meshreﬁnement driver. So far as I know it’s never been tested with any other driver.
 AHFinderDirect works ﬁne in single or multiprocessor Cactus runs.
 Obviously, your Cactus conﬁguration must include AHFinderDirect, and your ActiveThorns
parameter(s) must activate it.
 AHFinderDirect inherits from the other thorns (strictly speaking, implementations) listed in
table 1, so you’ll need them (or more precisely some thorns providing them) in your conﬁguration
and activated, too.
Implementation  Typically provided by Thorn 


Grid  CactusBase/CartGrid3d 
IO  CactusBase/IOUtil 
ADMBase  CactusEinstein/ADMBase 
StaticConformal  CactusEinstein/StaticConformal 
SpaceMask  CactusEinstein/SpaceMask 
SphericalSurface  AEIThorns/SphericalSurface 
Table 1: This table lists all the other implementations from which AHFinderDirect inherits, and
the thorns which typically provide these implementations.
 Grid::domain = "full", "bitant", "quadrant", and "octant" are supported. Alas, at present rotating
(or more precisely nonlocal) symmetry boundary conditions aren’t supported.
 The ADMBase::metric_type values "physical" and "static conformal" are supported; for the latter
you must have storage turned on for at least the conformal factor StaticConformal::psi. (The Cactus
3D grid functions for 1st and 2nd derivatives of psi aren’t used.)
 AHFinderDirect uses the CCTK_InterpGridArrays() Cactus global (multiprocessor grid array)
interpolator API; this is provided by PUGHInterp or CarpetInterp (so you must have the appropriate
one of these thorns compiled in and activated). CCTK_InterpGridArrays() in turn uses the new
CCTK_InterpLocalUniform() processorlocal interpolator API; AHFinderDirect uses various options in
this API which at present are only supported by thorn AEILocalInterp (so you must have this thorn
compiled in and activated)
 AHFinderDirect uses various Cactus reduction APIs to coordinate multiprocessor horizon ﬁnding, so
(even if you’re only going to run on a single processor) you must have a reduction thorn like
PUGHReduce or CarpetReduce compiled in and activated.
 At present only a few of AHFinderDirect’s parameters are steerable. (This is actually quite easy to ﬁx; I
just haven’t gotten around to it yet.)
 I think AHFinderDirect will “work” with checkpoint/restart, but I haven’t tested this yet. Here “work”
means the restart will be like starting a new run, in that AHFinderDirect will set the initial guess for
each horizon in a startofanewrun manner. Alas, it will also write a new BH_diagnostics ﬁle for
each horizon found, overwriting any existing BH_diagnostics ﬁles. This is a bug, which I
plan to ﬁx soon (the right behavior is/will be to append to the existing BH_diagnostics
ﬁle).
AHFinderDirect can pass information to the rest of Cactus in several ways; these are described in detail in
section 4.7.
3 Obtaining and Compiling AHFinderDirect
You should be able to obtain the source code for this thorn via the usual procedures for anonymous git checkout;
at present it lives in the EinsteinAnalysis arrangement.
This thorn is written primarily in C++, calling C and Fortran 77 numerical
libraries. In
theory the code should be quite portable to modern C++ compilers, but in practice I’ve had a
number of portability problems with various compilers. See the “Code Notes” and “Compiler Notes”
sections in the toplevel README ﬁle for details and lists of compilers currently known to be ok or
not.
By default AHFinderDirect doesn’t use any external libraries. However, if HAVE_DENSE_JACOBIAN__LAPACK is
deﬁned in src/include/config.h, then AHFinderDirect uses the LAPACK library for solving linear
equations. In this case you need to conﬁgure Cactus with LAPACK=yes. See the toplevel README and
README.library ﬁles for details on this.
4 AHFinderDirect Parameters
This thorn has lots of parameters, but most of them have reasonable default values which you probably won’t
need to change. Here I describe the parameters which you are likely to want to at least look at, and possibly set
explicitly.
Note that all of the “[$n$]”
parameters are Cactus array parameters, which you need to specify separately in your parameter ﬁle for each
apparent horizon. IMPORTANT: Apparent horizons are numbered starting at 1, not 0! The example
in section 8 should make this clear.
4.1 Overall Parameters

find_every
This is an integer parameter specifying how often AHFinderDirect should try to ﬁnd apparent
horizons: If $$find_every = 0,
AHFinderDirect is a noop. If $$find_every≠0,
AHFinderDirect tries to ﬁnd apparent horizons each find_every time steps.
The default value is 1, i.e. AHFinderDirect tries to ﬁnd apparent horizons at every time step.

N_horizons
How many apparent horizons do you want to ﬁnd in each slice? Typical values are 1 (the default), 2,
or 3.
This thorn numbers the apparent horizons from 1 to N_horizons inclusive. There are a
number of other parameters (described below) which you need to set for of these each apparent
horizons.
Note that N_horizons sets the number of apparent horizons you want to ﬁnd in the Cactus 3D
numerical grid, not in the whole spacetime. For example, if you are simulating (say) Misner data
with Grid::domain = "bitant", with the two throats at (say) roughly $z=\pm 1$,
then you should set N_horizons = 1 to ﬁnd those two apparent horizons, since you’re only ﬁnding
one apparent horizon within the numerical grid. If you also want to search for a common apparent
horizon surrounding both black holes, then you should set N_horizons = 2, since you’re ﬁnding at
most 2 apparent horizons within the numerical grid.

verbose_level
This controls how verbose this thorn is in printing informational (nonerror) messages describing what it’s
doing. In order from tersest to most verbose, the allowable values are

"no output"
Don’t print anything.

"physics highlights"
Print only a single line each time AHFinderDirect runs, giving which horizons were found.

"physics details"
Print two lines for each horizon found, giving the horizon area, centroid position, and
irreducible mass. This is the default.

"algorithm highlights"
Also print a single line for each Newton iteration giving the 2norm and $\infty $norm
of the $\Theta \left(h\right)$
function deﬁned by equation $\left(1\right)$.

"algorithm details"
Print lots of detailed messages tracing what the code is doing.

"algorithm debug"
Print even more detailed messages tracing what the code is doing, mainly useful for debugging
purposes.
4.2 Choosing the Local Coordinate Origin for each Apparent Horizon
For each apparent horizon you want to ﬁnd, you need to specify the Cactus
$\left(x,y,z\right)$
coordinates of a local coordinate system origin. As described in section 2, each apparent horizon must be a
Strahlkörper with respect to its local coordinate system origin.
You specify the local coordinate system origin for each horizon with the (Cactus array) parameters

origin_x[$n$] 
origin_y[$n$] 
origin_z[$n$] 

These all default to 0.0. In practice, you should set these parameters to be somewhere reasonably
close to your best guess for the center of each apparent horizon. These aren’t too critical: being oﬀ
by 1/4 of the horizon radius is no problem, and – assuming the algorithm still converges – even 1/2
of the horizon radius only slows the convergence by an extra iteration or two. But poor values of
these parameters do make the algorithm more likely to fail to converge.
At present the local coordinate origin is ﬁxed once you set it; there’s no provision
for it to move to track a moving black hole. I hope to add such a provision
soon.
4.3 Specifying the Initial Guess
AHFinderDirect requires an initial guess for the apparent horizon’s coordinate position and shape (that is, for
the $h\left(angle\right)$
function deﬁned in section 2), for each apparent horizon you want to ﬁnd. Unlike some other apparent
horizon ﬁnders (eg. the curvature ﬂow method in AHFinder), for AHFinderDirect there’s no
restriction on whether the initial guess is inside, outside, or crossing the actual apparent horizon: the
only important thing is that it should be “close”. Just how close the initial guess needs to be for
AHFinderDirect to ﬁnd the (a) apparent horizon depends on the slice and the coordinates, but as a general
rule of thumb any initial guess that’s within 1/3 to 1/2 of the mean horizon radius will probably
work.
The “initial guess” speciﬁcation is used the ﬁrst time we try to ﬁnd any given apparent horizon, and also any
succeeding time when the most recent attempt to ﬁnd this apparent horizon failed. If we succeed in ﬁnding a
given apparent horizon, than that apparent horizon position is automatically reused as the initial guess the next
time we try to ﬁnd the same apparent horizon; in this case the explicit “initial guess” speciﬁcation is
ignored.
There are a number of parameters for specifying the initial guess:

initial_guess_method[$n$]
This sets what type of the initial guess is used for each apparent horizon
position. There are several possibilities, most with their own sets of
subparameters: ${}^{,}$

"read from file"
This reads the initialguess $h\left(angle\right)$
function from a data ﬁle. The ﬁle format is currently hardwired to be that written with
file_format = "ASCII (gnuplot)" (see below). The subparameter initial_guess__read_from_named_file__file_name
speciﬁes the ﬁle name.

"Kerr/Kerr"
This sets the initial guess to the analyticallyknown apparent horizon position
in a Kerr spacetime in Kerr coordinates. This is a coordinate sphere of radius
$\left(1+\sqrt{1{a}^{2}}\right)m$.
There are subparameters

initial_guess__Kerr_Kerr__x_posn[$n$] 
initial_guess__Kerr_Kerr__y_posn[$n$] 
initial_guess__Kerr_Kerr__z_posn[$n$] 

to set the position of the Kerr black hole (note this position is in global Cactus coordinates,
not relative to the local coordinate origin), and

initial_guess__Kerr_Kerr__mass[$n$] 
initial_guess__Kerr_Kerr__spin[$n$] 

to set its mass and spin.

"Kerr/KerrSchild"
This sets the initial guess to the analyticallyknown apparent horizon position in a Kerr spacetime in
KerrSchild coordinates. This is a coordinate ellipsoid with radia (semimajor axes)
$$\begin{array}{ccc}{r}_{z}\hfill & \hfill =\hfill & \left(1+\sqrt{1{a}^{2}}\right)m\hfill \\ {r}_{x}={r}_{y}\hfill & \hfill =\hfill & {r}_{z}\sqrt{1+{\left(\frac{am}{{r}_{z}}\right)}^{2}}\hfill & \hfill =\hfill & \sqrt{\frac{2{r}_{z}}{m}}\phantom{\rule{0.3em}{0ex}}m\hfill \end{array}$$  (2) 

initial_guess__Kerr_KerrSchild__x_posn[$n$] 
initial_guess__Kerr_KerrSchild__y_posn[$n$] 
initial_guess__Kerr_KerrSchild__z_posn[$n$] 

(note this position is in global Cactus coordinates, not relative to the local coordinate
origin), and

initial_guess__Kerr_KerrSchild__mass[$n$] 
initial_guess__Kerr_KerrSchild__spin[$n$] 

to set its mass and spin.

"coordinate sphere"
This sets the initial guess to a coordinate sphere; there are subparameters

initial_guess__coord_sphere__x_center[$n$] 
initial_guess__coord_sphere__y_center[$n$] 
initial_guess__coord_sphere__z_center[$n$] 

(note this position is in global Cactus coordinates, not relative to the local coordinate
origin), and

initial_guess__coord_sphere__radius[$n$]

to set the radius.

"coordinate ellipsoid"
This sets the initial guess to a coordinate ellipsoid; there are subparameters

initial_guess__coord_ellipsoid__x_center[$n$] 
initial_guess__coord_ellipsoid__y_center[$n$] 
initial_guess__coord_ellipsoid__z_center[$n$] 

(note this position is in global Cactus coordinates, not relative to the local coordinate
origin), and

initial_guess__coord_ellipsoid__x_radius[$n$] 
initial_guess__coord_ellipsoid__y_radius[$n$] 
initial_guess__coord_ellipsoid__z_radius[$n$] 

to set the radia (semimajor axes).
4.4 I/O Parameters for the Apparent Horizon Shape(s)
The main output of this thorn is the computed horizon shape function
$h\left(angle\right)$, and
correspondingly the $\left(x,y,z\right)$
coordinate positions of the apparenthorizonsurface (angular) grid points. There are several parameters
controlling if, how often, and how these should be written to data ﬁles:

output_h_every
As described in section 4.1, AHFinderDirect will try to ﬁnd the apparent horizon(s) every
find_every time steps. However, you can control how often (if at all) the apparent horizon shape(s)
are written to data ﬁles: this is only done if output_h_every is nonzero, and the Cactus time step
number cctk_iteration is an integral multiple of this parameter (output_h_every).

file_format
This speciﬁes the ﬁle format for horizonshape (and other angulargridfunction) data ﬁles. Unfortunately,
at the moment only a single format is implemented,

"ASCII (gnuplot)"
This is a simple ASCII format designed for easy plotting with gnuplot:
 AHFinderDirect writes a separate data ﬁle for each Cactus time step and for each
apparent horizon found. By default these are all written in to the IOUtil::out_dir
directory; see h_directory (below) to change this.
 The time step number and the apparent horizon number are both encoded in the ﬁle
name; the actual ﬁle name is given by a printf() format "%s/%s.t%d.ah%d.%s",
where
 the ﬁrst %s is the directory set by the IO::out_dir and/or h_directory parameters
(see below)
 the second %s is the base ﬁle name set by the h_base_file_name parameter (see
below)
 the ﬁrst %d is the Cactus time step number
 the second %d is the apparent horizon number
 the third %s is the ﬁle name extension, set by the
ASCII_gnuplot_file_name_extension parameter; this defaults to "gp"
 Comment lines begin with #.
 Patches are separated by 2 blank lines; rows of apparenthorizon points within a patch are
separated by single blank lines.
 Each apparenthorizonsurface point is described by a single line, containing the
whitespaceseparated ﬁelds
 Two “unwrapped” angular coordinates in degrees, representing a point on ${S}^{2}$.
 The $h$
value, giving the radius of the apparent horizon surface at that angle.
 The corresponding Cactus $\left(x,y,z\right)$
coordinates.
See section 6 for a discussion of visualization for these and other AHFinderDirect output
ﬁles.

h_directory
This speciﬁes the directory in which the $h$
data ﬁles are to be written. If it doesn’t already exist, this directory is created before writing the data ﬁles.
This parameter defaults to the value of the IO::out_dir parameter.

h_base_file_name
This speciﬁes the base ﬁle name for $h$
data ﬁles, as described above. This defaults to ”h”.

ASCII_gnuplot_file_name_extension
This speciﬁes the ﬁle name extension for $h$
data ﬁles, as described above. This defaults to ”gp”.
4.5 Parameters for the “BH_diagnostics” Files
As well as the apparent horizon shape ﬁles, this thorn can also write ﬁles giving time series of various
diagnostics. These are controlled by the following parameters:

output_BH_diagnostics
If this Boolean parameter is set to true, AHFinderDirect will write a “black hole diagnostics” ﬁle for
each distinct apparent horizon found (up to N_horizons ﬁles in all). Each such ﬁle contains a
time series of various diagnostics for all time steps when the corresponding apparent horizon
was found. The ﬁle format is again a simple ASCII format designed for easy plotting with
gnuplot:
 The apparent horizon number is encoded in the ﬁle name; the actual ﬁle name is given by a printf()
format "%s/%s.ah%d.%s", where
 the ﬁrst %s is the directory set by the IO::out_dir and/or BH_diagnostics_directory
parameters (see below)
 the second %s is the base ﬁle name set by the BH_diagnostics_base_file_name parameter
(see below); this defaults to "BH_diagnostics"
 the %d is the apparent horizon number
 the third %s is the ﬁle name extension, set by the BH_diagnostics_file_name_extension
parameter; this defaults to "gp"
 The ﬁle begins with a block of header comments (lines begining with #) describing the data
ﬁelds.
 Each time this apparent horizon is found, a single line is appended to the data ﬁle, containing various
whitespaceseparated ﬁelds. The the precise list of ﬁelds, see the header comments, or see the
function output() in src/driver/BH_diagnostics.cc. As of this writing the ﬁelds
are:
 the Cactus iteration number cctk_iteration
 the Cactus time coordinate cctk_time
 the Cactus $\left(x,y,z\right)$
coordinates of the apparent horizon centroid
 the minimum, maximum, and mean coordinate radia of the apparent horizon about the
local coordinate origin
 the minimum and maximum Cactus $x$,
$y$,
and $z$
coordinates of the apparent horizon surface
 the proper circumferences of the apparent horizon in the $xy$,
$xz$,
and $yz$
localcoordinate planes
 the $xz\u2215xy$
and $yz\u2215xy$
ratios of the proper circumferences
 the proper area of the apparent horizon, $A$
 the irreducible mass of the apparent horizon, $\sqrt{A\u221516\pi}$
 the areal radius of the apparent horizon, $\sqrt{A\u22154\pi}$

BH_diagnostics_directory
This speciﬁes the directory in which the black hole diagnostics data ﬁles are to be written. If it doesn’t
already exist, this directory is created before writing the data ﬁles. This parameter defaults to the value of
the IO::out_dir parameter.

BH_diagnostics_base_file_name
This speciﬁes the base ﬁle name for black hole diagnostics data ﬁles, as described above. This defaults to
"BH_diagnostics".

BH_diagnostics_file_name_extension
This speciﬁes the ﬁle name extension for black hole diagnostics data ﬁles, as described above. This defaults
to "gp".
4.6 (Excision) Mask Parameters
This thorn can optionally set a mask grid function (or functions) at each point of the Cactus grid,
to indicate where that point is with respect to the apparent horizon(s). This is usually used for
excision.

set_mask_for_all_horizons 
set_mask_for_individual_horizon[$n$] 

These Boolean parameters control whether AHFinderDirect should set a mask grid function(s):
If the (C) expression set_mask_for_all_horizons  set_mask_for_individual_horizon[i] is
true, then AHFinderDirect will set a mask grid function(s) for apparent horizon $i$.
All these parameters default to false (don’t set a mask).
If any of these parameters is set to true, you almost certainly also need to set the Boolean parameter
SpaceMask::use_mask to true to to turn on storage for the mask grid function(s)!
If it’s setting a mask(s), AHFinderDirect partitions the Cactus grid into 3 regions: an “inside”, a “buﬀer”,
and an “outside”. Typically the inner region is excised, but AHFinderDirect doesn’t itself do this: It just sets
the mask(s); you need to use some other thorn(s) to do the actual excision.
The 3 regions are deﬁned as follows: For a grid point a distance
$r\left[i\right]$ from
horizon $i$’s local coordinate
origin, with horizon $i$’s
radius in this same direction (again, measured from its local coordinate origin) being
${r}_{horizon}\left[i\right]$,
the regions are deﬁned by
$$$$ 
r[i] ≤ rextr[i]  for some i     ⇒  inner

r[i] > rr[i]  for all i  and  r[i] ≤ rr[i]  for some i  ⇒  buﬀer

r[i] > rr[i]  for all i     ⇒  outer

(3)
where
$$\begin{array}{ccc}{r}_{r}\hfill & \hfill =\hfill & \hfill \end{array}$$mask_radius_multiplier× rn + mask_radius_offset× Δxe
rr =rr + mask_buffer_thickness× Δxe
 (4) 
and where ${\Delta x}_{e}$
is the basegrid Cactus grid spacing (more precisely, the geometric mean of the base grid’s
$x$,
$y$, and
$z$ Cactus grid
spacings) .

mask_radius_multiplier 
mask_radius_offset 
mask_buffer_thickness 

These parameters are used to deﬁne the radia ${r}_{r}$
and ${r}_{r}$
in equation $\left(4\right)$
above.
Note that the sign convention here is that mask_radius_multiplier is multiplied by the horizon
radius, then mask_radius_offset (scaled by the Cactus grid spacing) is added. Thus for use
with excision (where the inner region – which will be excised – must be somewhat inside the
horizon), mask_radius_multiplier should be a positive real number slightly less than 1.0, and/or
mask_radius_offset a negative real number.
The default values for these parameters are
mask_radius_multiplier = 0.8
mask_radius_offset = 5.0
mask_buffer_thickness = 5.0

mask_is_noshrink
This Boolean parameter speciﬁes whether the inside and buﬀer regions should be prevented from ever
shrinking during a time evolution (this is the default), or whether they should be set independently from
one time step to the next (and thus allowed to either grow or shrink). More precisely, once a given grid
point has been classiﬁed as inside, buﬀer, or outside, AHFinderDirect executes the following
algorithm:

inside


buﬀer

 if  (mask_is_noshrink && ($$mask = inside_value))

 # this point was previously inside $\Rightarrow $ noop here 
 else  $$mask ←buffer_value 

outside

 if  (mask_is_noshrink && (($$mask = inside_value)  ($$mask = buffer_value))

 # this point was previously inside or buﬀer $\Rightarrow $ noop here 
 else  $$mask ←outside_value 

min_horizon_radius_points_for_mask
By default, AHFinderDirect sets the mask for each apparent horizon found. If we’re using mesh
reﬁnement, it’s possible for an apparent horizon to be found on a coarse grid, and the masked region to be
only a few grid points across on a ﬁne grid. This causes some other Cactus thorns (eg. LegoExcision) to
crash. :(
This parameter can be used to avoid this problem: For each apparent horizon that it ﬁnds,
AHFinderDirect only sets the mask on a given grid if
$${r}_{inner,min}\ge $$min_horizon_radius_points_for_mask∗ Δxcurrent,max
 (5) 
where ${r}_{inner,min}$ is the minimum
over all angles of ${r}_{r}$ as
deﬁned by equation $\left(4\right)$,
and ${\Delta x}_{current,max}$ is the maximum
of the Cactus $x$,
$y$,
and $z$ grid spacings on the
current Cactus grid.
If this condition isn’t satisﬁed, then AHFinderDirect skips setting the mask for this apparent horizon,
just as if this apparent horizon wasn’t found. (Note that all other processing for an apparent horizon being
found is still done, including writing output ﬁles, using the apparent horizon shape as the initial guess for
the next time step’s apparenthorizon ﬁnding, etc.; it’s only the mask processing for this horizon (at this
time step) that’s skipped.)
AHFinderDirect supports two types of mask grid functions; the following two Boolean parameters choose
which of them you want to set; you can set either or even both of these:

set_old_style_mask
This parameter (default true) speciﬁes an oldstyle excision mask, one stored in a CCTK_REAL
Cactus grid function. (The AHFinder apparent horizon ﬁnder uses this type of mask.)

set_new_style_mask
This parameter (default false) speciﬁes a newstyle excision mask, one stored in a speciﬁed bit
ﬁeld of a CCTK_INT Cactus grid function. The bit ﬁeld is speciﬁed by its name, as registered
with the SpaceMask thorn. We plan to eventually convert all Cactus excision (and other uses
of mask grid functions) to this scheme, but at the moment not much code supports it. Note that
AHFinderDirect doesn’t itself create/register any bit ﬁelds or state names with SpaceMask –
you must arrange for some other thorn(s) to do this.
For an oldstyle mask, the following parameters specify the mask grid function and how it should be
set:

old_style_mask_gridfn_name
This parameter speciﬁes the mask grid function’s name.

old_style_mask_inside_value 
old_style_mask_buffer_value 
old_style_mask_outside_value 

If an oldstyle mask is to be set in the corresponding regions, these parameters specify the values
to which it should be set. These are all CCTK_REAL values.
For an newstyle mask, the following parameters specify the mask grid function and how it should be
set:

new_style_mask_gridfn_name 
new_style_mask_bitfield_name 

These parameters specify the mask grid function’s name and the bitﬁeld name within it.

new_style_mask_inside_value 
new_style_mask_buffer_value 
new_style_mask_outside_value 

If an newstyle mask is to be set in the corresponding regions, these parameters specify the values to
which it should be set. These are all characterstring state names, as registered with the SpaceMask
thorn.
Note that AHFinderDirect doesn’t itself register any bitﬁelds or states with SpaceMask – you must arrange
for some other thorn(s) to do this before AHFinderDirect tries to ﬁnd the horizon(s).
If AHFinderDirect sets a mask or masks, this happens in the same schedule bin(s) as the horizon ﬁnding.
More precisely, AHFinderDirect creates two schedule groups for this purpose:
 The schedule group group_for_mask_stuff is scheduled to run just after the horizon is found.
 The schedule group group_where_mask_is_set is scheduled inside the schedule group
group_for_mask_stuff.
 The actual setting of the mask is scheduled inside the schedule group group_where_mask_is_set.
Thorn PreviousMask uses these schedule groups to keep a “previous” as well as a “current” mask. See that
thorn’s thorn guide for further details.
4.7 Communicating with Other Thorns
Besides the data ﬁles it writes, AHFinderDirect currently has three ways to communicate with other Cactus
thorns:
 AHFinderDirect can set a mask grid function(s) based on the (a) horizon’s shape; other thorns
may then use this for excision or other purposes. AHFinderDirect supports both the oldstyle
(CCTK_REAL) mask (compatible with AHFinder) or the newstyle (CCTK_INT) mask bitﬁelds
deﬁned by SpaceMask; you can even use both styles simultaneously.
 AHFinderDirect can announce a selected horizon’s centroid position to another thorn (typically
DriftCorrect, which uses this to adjust its corotating shift vector). This uses the new
functionaliasing version of DriftCorrect, not the old version which worked with an auxiliary thorn
AHFSetDCCentroid.
 AHFinderDirect provides a set of aliased functions which any other thorn(s) can call to ﬁnd out
the shape of a speciﬁed horizon.
 AHFinderDirect can store information about the horizon(s) it ﬁnds in the SphericalSurface
variables for other thorns to use.
AHFinderDirect’s mask features are described in section 4.6; the other communication mechanisms are
described in the following subsections.
4.7.1 Parameters for Announcing a Horizon Centroid to Other Thorns
This thorn can optionally announce the centroid of a speciﬁed apparent horizon to another thorn
(typically DriftCorrect) each time that apparent horizon is found. This is controlled by the following
parameter:

which_horizon_to_announce_centroid
This is an integer parameter which defaults to 0 (which means not to announce any centroid). If
it’s set to a nonzero integer, that speciﬁes the horizon number to have its centroid announced.
4.7.2 Aliased Functions to Provide HorizonShape Information
AHFinderDirect provides the following aliased functions to allow other thorns
to ﬁnd out about the horizons. Each function returns a status code which is
$\ge 0$ for ok,
or negative for an error.
#
# This function returns the local coordinate origin for a given horizon.
#
CCTK_INT FUNCTION HorizonLocalCoordinateOrigin \
(CCTK_INT IN horizon_number, \
CCTK_REAL OUT origin_x, CCTK_REAL OUT origin_y, CCTK_REAL OUT origin_z)
#
# The following function queries whether or not the specified horizon
# was found the most recent time AHFinderDirect searched for it.
# The return value is:
# 1 if the horizon was found
# 0 if the horizon was not found
# negative for an error
#
CCTK_INT FUNCTION HorizonWasFound(CCTK_INT IN horizon_number)
#
# The following function computes the horizon radius in the direction
# of each (x,y,z) point, or 1.0 if this horizon wasn’t found the most
# recent time AHFinderDirect searched for it. More precisely, for each
# (x,y,z), consider the ray from the local coordinate origin through
# (x,y,z). This function computes the Euclidean distance between the
# local coordinate origin and this ray’s intersection with the horizon,
# or 1.0 if this horizon wasn’t found the most recent time AHFinderDirect
# searched for it.
#
# If this function is to be used in a multiprocessor run on a horizon
# which was found on some other processor, the parameter
# AHFinderDirect::always_broadcast_horizon_shape
# must be set to true to get the correct answer.
#
CCTK_INT FUNCTION HorizonRadiusInDirection \
(CCTK_INT IN horizon_number, \
CCTK_INT IN N_points, \
CCTK_REAL IN ARRAY x, CCTK_REAL IN ARRAY y, CCTK_REAL IN ARRAY z, \
CCTK_REAL OUT ARRAY radius)
4.7.3 Storing HorizonShape Information in the SphericalSurface Variables
SphericalSurface (in the AEIThorns arrangement) deﬁnes a set of generic grid arrays which describe
“spherical surfaces”. AHFinderDirect can optionally store information about the horizons it ﬁnds in the
SphericalSurface variables. This is controlled by the following parameters:

which_surface_to_store_info[$n$]

This parameter should be set to the SphericalSurface surface number into which information on
a given AHFinderDirect horizon should be stored, or to 1 to skip storing the information. It
defaults to 1 for each horizon.
Note that SphericalSurface numbers surfaces starting from 0, whereas AHFinderDirect
numbers horizons starting from 1!
At present, if multiple AHFinderDirect horizons specify the same SphericalSurface surface,
the highestnumbered horizon will “win”, i.e. it will overwrite the data from any lowernumbered
horizons.
4.8 Other Parameters

run_at_CCTK_ANALYSIS 
run_at_CCTK_POSTSTEP 
run_at_CCTK_POSTINITIAL 
run_at_CCTK_POST_RECOVER_VARIABLES 

These parameters (which default to true, false, false, and true respectively) control which
schedule bins AHFinderDirect runs in. Historically, AHFinderDirect ran in CCTK_ANALYSIS,
and that’s still the default, but these parameters allow you to change this so it runs in CCTK_POSTSTEP
and/or CCTK_POSTINITIAL instead. (You can even run in all three bins if you want!)
In general we need to run at CCTK_POST_RECOVER_VARIABLES, since
 parameters may have been steered at recovery, so we may need to ﬁnd a new horizon or
horizons, and
 we need to set the mask again to make sure it’s correct right away (since our next regular
horizonﬁnding may not be until some time steps later)
Therefore the run_at_CCTK_POST_RECOVER_VARIABLES parameter should probably be left at its default
setting of true.

geometry_interpolator_name 
geometry_interpolator_pars 

These parameters control the (3D) “geometry interpolation” of the spacetime geometry
(${g}_{ij}$ and
${K}_{ij}$, or
their StaticConformal equivalents) to the apparent horizon position. The defaults are set to use a quadratic
Hermite interpolator. This works fairly well, but because of the interpolator molecule size you must use
$$driver::ghost_size ≥ 2.
If you want to get very high accuracy from AHFinderDirect, then you should use a cubic Hermite
geometry interpolator, by setting
AHFinderDirect::geometry_interpolator_pars = "
order=3
boundary_off_centering_tolerance={1.0e10 1.0e10 1.0e10 1.0e10 1.0e10 1.0e10}
boundary_extrapolation_tolerance={0.0 0.0 0.0 0.0 0.0 0.0}
"
Assuming perfectly accurate geometry variables in the 3D Cactus grid, this will make
AHFinderDirect (very) roughly an order of magnitude more accurate. However, the larger
molecule size will make it about a factor of 2–3 slower, and will also require that you set
$$driver::ghost_size ≥ 3.
The sample parameter ﬁle par/Kerrorder3.par shows an example of this.

N_zones_per_right_angle[$n$]
This parameter sets the angular resolution used to compute each patch. The units are the number of
angular grid zones per right angle. The default is 18, i.e. a 5 degree angular resolution. There’s no
problem with this parameter varying from one horizon to another, but for simplicity it should be
even.
For any horizon which is close to spherical about its local coordinate origin, you can lower this parameter
to make AHFinderDirect run faster (typical runtimes scale roughly as the cube of this parameter); 6 is
about the minimum reasonable value.
For any horizon which is highly nonspherical about its local coordinate origin, you can raise
this parameter to get better resolution; 30 should be enough for even a highly nonspherical
horizon.

max_allowable_horizon_radius[$n$]
This parameter gives the maximum meancoordinateradius which any given
trial surface may have in the course of trying to solve the apparent horizon
equation.
In particular, if any trial surface has a mean coordinate radius which exceeds this parameter,
AHFinderDirect gives up and deems this apparent horizon to be “not found”.
This parameter defaults to $1{0}^{10}$
(eﬀectively $+\infty $)
for each apparent horizon. You can set it to a smaller value to make AHFinderDirect a bit more
eﬃcient, or (probably more important in practice) to stop AHFinderDirect from iterating oﬀ the edge of
the grid if this causes problems with interpolation or boundary conditions.

max_allowable_Theta
This parameter gives the maximum $\parallel \Theta {\parallel}_{\infty}$
which any given trial surface may have in the course of trying to solve the apparent horizon equation. In particular, if any
trial surface has $\parallel \Theta {\parallel}_{\infty}$
exceeding this parameter, AHFinderDirect gives up and deems this apparent horizon to be “not found”. This parameter
defaults to $1{0}^{10}$
(eﬀectively $+\infty $)
for each apparent horizon.

surface_expansion[$n$]

This parameter (which defaults to 0.0) sets the expansion of each surface. With the default setting this thorn
solves $\left(1\right)$
to ﬁnd apparent horizons. With other settings of this parameter this thorn can be used to ﬁnd “surfaces of
constant expansion”; these may be useful for excision, wave extraction, or other purposes
([?]).
To help in choosing the value(s) of the
surface_expansion[$n$]
parameter, ﬁgure 1 (from [?]), shows the expansion of
$r=constant$
surfaces in an EddingtonFinkelsteon slice of the unitmass Schwarzschild spacetime.
5 Monitoring AHFinderDirect’s Status
There are two primary ways of monitoring what AHFinderDirect is doing during a Cactus run: the
BH_diagnostics ﬁles and the CCTK_INFO messages written to the Cactus standard output:
The BH_diagnostics ﬁles are described in detail in section 4.5. These ﬁles are written and “ﬂushed” at each
time step, so they’re always uptodate.
During the apparenthorizon–ﬁnding process, AHFinderDirect writes various CCTK_INFO messages describing
the convergence of the iterative solution of the apparent horizon equation 1 on each processor. In particular, if
verbose_level is set to "algorithm highlights" or a more verbose setting (cf. section 4.1), then
AHFinderDirect writes CCTK_INFO messages like these:
INFO (AHFinderDirect): proc 0/horizon 1:it 1 r_grid=0.595 Theta=1.1e01
INFO (AHFinderDirect): proc 0/horizon 1:it 2 r_grid=0.614 Theta=7.2e02
INFO (AHFinderDirect): proc 0/horizon 1:it 3 r_grid=0.632 Theta=2.9e02
INFO (AHFinderDirect): proc 0/horizon 1:it 4 r_grid=0.642 Theta=9.9e04
INFO (AHFinderDirect): proc 0/horizon 1:it 5 r_grid=0.642 Theta=7.9e07
INFO (AHFinderDirect): proc 0/horizon 1:it 6 r_grid=0.642 Theta=7.2e13
INFO (AHFinderDirect): AH 1/2: r=0.660716 at (0.000000,0.000000,1.127434)
INFO (AHFinderDirect): AH 1/2: area=338.0473838 m_irreducible=2.59330658
INFO (AHFinderDirect): writing h to "misner.h.t0.ah1.gp"
Here r_grid is a rough estimate of the mean radius of the trial surface at each iteration, and Theta is the
inﬁnitynorm of $\Theta $,
the left hand side of the apparent horizon equation 1 over the surface. Once the apparent
horizon has been found (Theta is suﬃciently small), then AHFinderDirect prints its mean
radius,
centroid position, area, and irreducible mass.
6 Visualization
There are several ways to visualize AHFinderDirect’s output:
The simplest is to plot various quantities from the BH_diagnostics ﬁles (described in detail in section 4.5). For
example, using gnuplot (http://www.gnuplot.info), you can plot a graph of the surface area of horizon #4 as
a function of coordinate time, with the command
plot ’BH_diagnostics.ah4.gp’ using 2:26 with points
ygraph (http://www.aei.mpg.de/~pollney/ygraph/) may also be able to directly plot the BH_diagnostics
ﬁles.
Given a horizonshape data ﬁle h.t105.h4.gp, the gnuplot command
splot ’h.t105.h4.gp’ with lines
will plot the $h\left(angle\right)$
function, with the $x$
and $y$ axes
of the plot being the two “unwrapped” angular coordinates on
${S}^{2}$, in degrees,
and the $z$ axis
being $h\left(angle\right)$.
However, in practice this usually isn’t very informative. Instead, you probably want the gnuplot
command
splot ’h.t105.h4.gp’ using 4:5:6 with lines
which will plot the 3D shape of the apparent horizon surface.
The src/misc directory of the AHFinderDirect source code contains several perl scripts which are
useful in visualizing AHFinderDirect output. In particular, the script select.plane selects a
particular 2D plane. If you put this in your Unix path, it can be used with a gnuplot command
like
splot ’<select.plane xy <h.t105.h4.gp’ using 4:5 with lines
to show the shape of an apparent horizon in the xy
plane.
Another Visualization option is OpenDX (http://www.opendx.org/). Thomas Radke’s has written some
OpenDX macros ImportAHFinderDirectGnuplot.net and ImportAHFinderDirectGnuplotPatch.net
to import AHFinderDirect horizonshape data ﬁles. These macros use a set of “control ﬁles”
named *.dx, one per horizon, which AHFinderDirect (by default) writes into the same directory
as the main horizon–shape output ﬁles. You can get these macros by anonymous CVS with the
command
cvs d :pserver:cvs_anon@cvs.aei.mpg.de:/numrelcvs \
checkout AEIPhysics/Visualization/OpenDX
Another Visualization option is to produce ﬁles in the standard Xdmf [?] ﬁle format which can be visualized for
example using VisIt [?]. Frank Löﬄer has written a python script AH2xdmf.py, which is included in this
thorn’s code repository, which generates such ﬁles from AHFinderDirect’s output from horizon shape ﬁles
h.t%d.ah%d.gp\. Please see its help text for details.
7 Accuracy
The apparent horizon positions are typically computed very accurately; tests on Kerr spacetimes give typical errors
of $1{0}^{4}m$ to
$1{0}^{5}m$.
The various diagnostics printed to standard output and written to the black hole diagnostics ﬁle(s), are typically
computed to accuracies on the order of a part per million or so.
Note, however, that the irreducible mass ${m}_{e}$
may diﬀer considerably from the black hole’s local mass or its contribution to the
slice’s ADM mass. For example, for Kerr spacetime in KerrSchild coordinates,
${m}_{e}\u2215{m}_{M}=0.949$,
$0.894$, and
$0.723$ for spin
parameters $a\equiv J\u2215{m}^{2}=0.6$,
$0.8$, and
$0.999$,
respectively. It would be better to (also) use the “isolated horizons” formalism of [?]; at some point this thorn
may be enhanced to do this.
8 Examples
There are a few example parameter ﬁles in the par/ directory, including Kerr initial data, Misner initial data,
and Misner timeevolution tests. The Kerrtiny.par parameter ﬁle is close to a minimal AHFinderDirect
example:
# This parameter file sets up Kerr/KerrSchild initial data, then
# finds the apparent horizon in it. The local coordinate system origin
# and the initial guess are both deliberately decentered with respect
# to the black hole, to make this a nontrivial test for the apparent
# horizon finder.
#
# This parameter file is "tiny" in the sense that it sets only a
# small number of AHFinderDirect parameters.
# flesh
cactus::cctk_itlast = 0
ActiveThorns = "PUGH"
driver::ghost_size = 2
driver::global_nx = 31
driver::global_ny = 31
driver::global_nz = 19
ActiveThorns = "CoordBase CartGrid3D"
grid::domain = "bitant"
grid::avoid_origin = false
grid::type = "byspacing"
grid::dxyz = 0.2
ActiveThorns = "ADMBase ADMCoupling StaticConformal Spacemask CoordGauge Exact"
ADMBase::initial_lapse = "exact"
ADMBase::initial_shift = "exact"
ADMBase::initial_data = "exact"
ADMBase::lapse_evolution_method = "static"
ADMBase::shift_evolution_method = "static"
ADMBase::metric_type = "physical"
Exact::exact_model = "Kerr/KerrSchild"
Exact::Kerr_KerrSchild__mass = 1.0
Exact::Kerr_KerrSchild__spin = 0.6
########################################
ActiveThorns = "IOUtil"
IOUtil::parfile_write = "no"
########################################
ActiveThorns = "SphericalSurface"
ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
AHFinderDirect::h_base_file_name = "Kerrtiny.h"
AHFinderDirect::N_horizons = 1
AHFinderDirect::origin_x[1] = 0.5
AHFinderDirect::origin_y[1] = 0.7
AHFinderDirect::origin_z[1] = 0.0
AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
AHFinderDirect::initial_guess__coord_sphere__x_center[1] = 0.2
AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.3
AHFinderDirect::initial_guess__coord_sphere__z_center[1] = 0.0
AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
9 Surfaces of Constant Expansion
Surfaces of Constant Expansion (CE surfaces) are introduced in [?] as a generalisation of apparent horizons
(AH). On an AH surface, the expansion is zero everywhere. On a CE surfaces, the expansion is still everywhere
the same, but it need not be zero. CE surfaces are also a generalisation of Constant Mean Curvature surfaces
(CMC surfaces); both are identical when the extrinsic curvature vanishes. As described in [?], it is likely that CE
surfaces foliate the spacelike hypersurface outside of some interior region. This interior region is inside the
common apparent horizon, if it exists.
CE surfaces can give some insight into the spacetime, because they can be used to analyse the part of the
spacelike hypersurface “between the horizons and inﬁnity”. Most notably, they can be used to look at the region
where a common horizon is about to (or believed to) form. Similarly, one can use them for collapsing stars where
an apparent horizon has not yet formed.
10 Pretracking
Apparent horizon pretracking is introduced in [?]. This is an application of CE surfaces. Even when there is no
common horizon, there are still common CE surfaces surrounding multiple black holes. Pretracking consists of
tracking in time the smallest common CE surface that can be found. It is reasonable to believe that this surface
will evolve into the common horizon at the time where this common horizon begins to exist. The expansion of
this smallest CE surface is also an indication of how close the spacelike hypersurface is to having a common
apparent horizon.
11 How AHFinderDirect Works
AHFinderDirect uses the apparent horizon (henceforth “horizon”) ﬁnding algorithm of [?], modiﬁed slightly to
work with ${g}_{ij}$
and ${K}_{ij}$ on a
Cartesian ($xyz$)
grid. The algorithm is described in detail in [?].
11.1 General Description of the Algorithm
As described above, I parameterizes the horizon shape by
$r=h\left(angle\right)$ for some singlevalue function
$h:{S}^{2}\to {\Re}^{+}$. The apparent horizon
equation $\left(1\right)$ then becomes
a 2D elliptic PDE on ${S}^{2}$
for the function $h$.
I ﬁnite diﬀerence this in angle to obtain a system of simultaneous nonlinear algebraic equations for
$h$ at the
angular grid points, and solve this system of equations by a global Newton’s method (or a variant with improved
convergence).
Computationally, this algorithm has 3 main parts:
11.2 The Multipatch System
Perhaps the most unusual feature of AHFinderDirect is the “multipatch” system used to cover
${S}^{2}$ without
coordinate singularities. In general there are 6 patches, one each covering a neighborhood of the
$\pm z$,
$\pm x$, and
$\pm y$ axes, but this may
be reduced in the presence of suitable symmetries. For example, ﬁgure 2 on page 42 shows a system of 3 patches covering
the $+xyz$ octant
of ${S}^{2}$.
This would be suitable for ﬁnding an apparent horizon with mirror symmetry about the (local)
$z=0$ plane,
and either 90 degree periodic rotation symmetry about the (local)
$z$ axis, or mirror symmetry
about each of the (local) $x$ and
$y$ axes.
To allow easy angular ﬁnite diﬀerencing within the patch system, each patch is extended beyond its nominal extent by a
“ghost zone”
(2 grid points wide in ﬁgure 2). Angular grid function values in the ghost zone can be obtained by interpatch
interpolation
or by applying symmetry operations. Once this is done, then angular ﬁnite diﬀerencing within the nominal
extent of each patch can proceed normally, ignoring the patch boundaries. AHFinderDirect can be
conﬁgured at compile time to use either 2nd order or 4th order angular ﬁnite diﬀerencing (3 point or
5 point angular molecules); the default is 4th order (5 point). This is conﬁgured at compile time in
src/include/config.h.
By default AHFinderDirect will automagically choose a patch system type for each apparent
horizon searched for, based on the local coordinate origin and the symmetries implicit in the Cactus
grid type. This generally works well, but if desired you can instead manually specify the patch
system type, the angular resolution, the width of the ghost zones, etc. See the param.ccl ﬁle for
details.
11.3 Other Software Used
AHFinderDirect’s src/sparsematrix/ directory contains various sparsematrix libraries, which have their
own copyrights and licensing terms:
The src/sparsematrix/umfpack/ directory contains a subset of the ﬁles in UMFPACK version 4.0
(11.Apr.2002). This code is copyright (©) 2002 by Timothy A. Davis, and is subect to the UMFPACK
License:
Your use or distribution of UMFPACK or any modified version of
UMFPACK implies that you agree to this License.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use or copy this program, provided
that the Copyright, this License, and the Availability of the original
version is retained on all copies. User documentation of any code that
uses UMFPACK or any modified version of UMFPACK code must cite the
Copyright, this License, the Availability note, and "Used by permission."
Permission to modify the code and to distribute modified code is granted,
provided the Copyright, this License, and the Availability note are
retained, and a notice that the code was modified is included. This
software was developed with support from the National Science Foundation,
and is provided to you free of charge.
12 Other Related Thorns
If you’re interested in AHFinderDirect, you might also be interested in some other related thorns:

EHFinder
 (in the AEIDevelopment arrangement) was written by Peter Diener, and ﬁnds the event
horizon(s) in a numerically computed spacetime. It’s described in detail in the paper [?].

AHFinder
 (in the CactusEinstein arrangement) was written by Miguel Alcubierre, and includes two
diﬀerent algorithms for ﬁnding apparent horizons, a minimization method and a “fast ﬂow” method
based on [?]. Unfortunately, both methods are very slow in practice.

TGRapparentHorizon2D
 (in the TAT arrangement) was written by Erik Schnetter, and is another
apparent horizon ﬁnder. It uses methods very similar to this thorn, and (like this thorn) is very
fast and accurate. However, it’s no longer under active development. It’s described in detail in the
papers [?] and [?].