Cactus Code Thorn AHFinder
Author(s) : Miguel Alcubierre
Maintainer(s): Cactus team
Licence : LGPL
--------------------------------------------------------------------------
1. Purpose
This thorn looks for apparent horizons (AH) in 3D. An AH is defined
as a surface where the expansion of outgoing null geodesics is zero.
NOTE: THIS FILE IS SOMEWHAT OUT OF DATE. I WILL UPDATE IT WHEN
I HAVE MORE TIME. PLEASE LOOK AT THE FILE `MinimumAHF_param.h'
FOR A MORE UP TO DATA DESCRIPTION OF THE PARAMETERS.
BRIEF DESCRIPTION
-----------------
This thorn looks for apparent horizons (AH) in 3D using two
different algorithms. In the default mode, it uses a minimization
algorithm. It can also use a flow algorithm if one sets the
parameter "ahf_flow" equal to "yes":
ahf_flow = "yes"
The minimization algorithm is VERY slow if one has a lot of terms
in the spherical harmonics expansion (see below), but it is quite
fast for a small number of terms.
The thorn is activated by adding the following line to the parameter
file:
ahfinder::ahf_active = "yes"
Also, one needs to tell the routine how often we want to
look for AH by using the parameter `mahf_findevery' (its
default value is 1).
The routine defines the grid function `mahfgrid' as:
ahfgrid = r - f(theta,phi)
So that the surface being considered corresponds to the zero level of
`mahfgrid'. The function f(theta,phi) is expanded in spherical
harmonics Y(l,m). Notice, however, that I DO NOT use the standard
normalization factor:
sqrt( (2l + 1) (l - m)! / 4 pi (l + m)! )
I use instead:
sqrt( 2l + 1 ) m = 0
sqrt( 2 (2l + 1) (l - m)! / (l + m)! ) m != 0
Notice first that I leave out the 4*pi. This is because I want Y(0,0)
to be the real radius of a sphere. Also, for M non-zero I have an
extra sqrt(2). I need this because I use a real basis with sines and
cosines instead of complex exponentials.
The minimization algorithm has three main modes of operation:
* Default mode: In the default mode, the routine tries to
find a minimum of the integral of the square of the expansion.
At an AH, this integral should be zero, and for any other
surface it will be positive. One must remember that the
routine can easily end up in a local minimum that is not
a horizon. To try to prevent this, the routine first
looks through a subset of the parameter space to try
to find a good guess for the global minimum. This, of
course can fail.
* Area minimization mode: This mode is activated by
setting the parameter:
mahf_minarea = "yes"
In this mode, the routine looks for a local minimum of the surface
area. For time-symmetric data (and only then), such a minimum will
correspond to an AH.
* Trapped surface finding mode: This mode is activated by
setting the parameter:
mahf_trapped_surface = "yes"
In this mode, the routine looks for a local minimum of the integral of
the square of the quantity (expansion + trapped_surface_delta). If
trapped_surface_delta is positive, then the routine will try to find a
surface whose expansion is -(trapped_surface_delta) everywhere on the
surface. One can monitor the number of interpolated points that have
a negative expansion by setting veryverbose = "yes", which indicates
whether or not a trapped surface has truly been found.
OUTPUT
------
The routine controls its own output. It produces output of 2 grid
functions (mahfgrid, mahf_exp) on its own (so please don't add them to
the Cactus output). The reason for this is that the finder might be
called at times that have nothing to do with when Cactus wants to do
output.
The grid function `mahfgrid' was described above. The grid function
`mahf_exp' is the expansion of outgoing photons on the level sets of
`mahfgrid'.
The finder also produces a series of data files:
1) "mahf_coeff.alm"
This file contains the coefficients of the spherical
harmonics expansion of the surface found. Notice
that the finder almost always finds a surface, but this
might not be a horizon, so be careful.
The format in which these coefficients are saved is the following:
* Lines beginning with # are comments. These include the time
at which the finder was called, and a line stating if the
surface found is a horizon or not.
* Empty lines indicate different calls to the finder (used when
the finder is called many times during an evolution).
* The coefficients are saved using a loop of the form:
do l=1,lmax
do m=-l,l
write(*,*) a(l,m),l,m
end do
end do
Notice that because our functions are real, I use an expansion
in sines and cosines and I really only consider positive values
of m. On output my convention is that coefficients with positive
m represent the terms:
P_(l,m) cos(m theta)
and coefficients with negative m represent the terms:
P_(l,m) sin(m theta)
2) "mahf.gauss"
This file contains a map of the gaussian curvature on the surface.
The format of this file might still change in the near future.
I need to talk to Werner about this.
* Lines beginning with # are comments.
* The data is written is a loop:
do i=0,ntheta-1
do j=0,nphi-1
write gaussian(i,j)
end do
end do
theta and phi are subdivided uniformly
(according to grid type) in the intervals:
octant: theta=[0,pi/2] phi=[0,pi/2]
quadrant: theta=[0,pi] phi=[0,pi/2]
full: theta=[0,pi] phi=[0,2 pi]
3) "mahf_area.tl"
Area of the surface for each time the finder was called.
Again, remember that the surface might not be a horizon.
4) "mahf_mass.tl"
Mass of the surface for each time the finder was called.
This is defined as:
M = sqrt(A/(16 pi))
Again, remember that the surface might not be a horizon.
5) "mahf_circ_eq.tl"
Equatorial circumference of the surface.
Again, remember that the surface might not be a horizon.
6) "mahf_meri_p1.tl"
Length of meridian of surface for phi=0.
Again, remember that the surface might not be a horizon.
7) "mahf_meri_p2.tl"
Length of meridian of surface surface for phi=pi/2.
Again, remember that the surface might not be a horizon.
8) "mahf_logfile"
Log file for the last time the horizon was called. The
reason why only the last call is here is that this file
is very long.
PARAMETERS
----------
Other parameters for the thorn are (see also file MinimumAHF_param.h):
STRINGS:
mahf_logfile = [yes,no] Write a log file. If yes, the
thorn produces the file "mahf_logfile".
(default = "no")
mahf_verbose = [yes,no] Write messages to screen.
(default = "yes")
mahf_veryverbose = [yes,no] Write lots of messages to the screen,
essentially the whole log file.
(default = "no")
mahf_guessverbose = [yes,no] Write a little info for each trial point
calculated in the guess process.
(default = "no")
mahf_2Doutput = [yes,no] 2D output of grid functions?
(default = "no")
mahf_3Doutput = [yes,no] 3D output of grid functions?
(default = "no")
mahf_mask = [yes,no] Use mask for definite horizons?
(default = "no")
mahf_weakmask = [yes,no] Use mask for possible horizons?
(default = "no")
mahf_phi = [yes,no] Expand in phi. It is useful to switch
this off if one knows in advance that
the horizon is really axisymmetric.
(default = "yes")
mahf_offset = [yes,no] Is the center offset from origin?
(default = "no")
mahf_wander = [yes,no] Do we allow the center to wander?
(default = "no")
mahf_refx = [yes,no] Do we have reflection symmetry x -> -x?
(default = "no")
mahf_refy = [yes,no] Do we have reflection symmetry y -> -y?
(default = "no")
mahf_refz = [yes,no] Do we have reflection symmetry z -> -z?
(default = "no")
mahf_octant = [yes,no,high] Does the surface have octant symmetry?
The "high" option forces also reflection
symmetry across the x-y diagonal.
(default = "no")
mahf_areamap = [yes,no] Construct an area map?.
(default = "no")
mahf_sloppyguess = [yes,no] This considers only spheres for the
initial guess. It is much faster.
(default = "no")
mahf_guessold = [yes,no] Use horizon found in previous call as
initial guess? Only relevant for
evolutions.
(default = "no")
mahf_inner = [yes,no] Tries to look for inner horizon instead.
(default = "no")
mahf_guess_absmin = [yes,no] Use the absolute minimum of the surface
integral of the square of the expansion
as the place to start the minimization
process. This is useful if you have
a very dynamical spacetime, with possibly
many local minima and you don't have any
idea where (or even if) you have an AH.
(default = "no")
mahf_manual_guess = [yes,no] Use a manually specified guess, via. the
parameters mahf_l0_guess, mahf_l2_guess,
etc. (works only in octant mode right now)
(default = "no")
INTEGERS:
mahf_findafter After how many timesteps start looking
for horizons (default 0).
mahf_findevery How often to look for horizons (iterations).
(default 1)
mahf_maxiter Maximum number of iterations of
Powell's minimization algorithm.
(default 10)
mahf_flowiter Maximum number of iterations for flow
(default 200)
mahf_lmax (<20) Number of terms in theta expansion
(default 2)
mahf_ntheta Number of subdivisions in theta for
surface integrals (default 200).
mahf_nphi Number of subdivisions in phi for
surface integrals (default 200).
mahf_nn0 Number of subdivisions of c0(0) for
initial guess (default 10)
mahf_nn2 Number of subdivisions of c0(2) for
initial guess (default 10)
REALS:
mahf_findaftertime After what time start looking for horizons.
A non-zero value of this parameter overrides
the paramter `mahf_findafter' above.
(default 0.0)
mahf_r0 Radius of initial sphere (0 forces largest
possible sphere)
(default 0.0)
mahf_xc x coordinate of center of expansion
(default 0.0)
mahf_yc y coordinate of center of expansion
(default 0.0)
mahf_zc z coordinate of center of expansion
(default 0.0)
mahf_tol Fractional tolerance for minimization
algorithm. If we decrease by less than
this in one iteration we are done.
(default = 0.1)
trapped_surface_delta In 'mahf_trapped_surface' mode, this
determines the expansion of the surface
that one is looking for. Notice that a
positive value of 'mahf_trapped_surface'
will cause the finder to look for a
surface with expansion equal to MINUS
this value everywhere.
(default = 0.0)
mahf_flowa Alpha parameter for flow (default 0.01).
mahf_flowb Beta parameter for flow (default 0.01).
mahf_flowh Weight of H flow (default 0.0).
mahf_flowc Weight of C flow (default 1.0).
mahf_flown Weight of N flow (default 0.0).
mahf_flowtol Tolerance for flow (default 0.0001).
mahf_maskshrink Shrink factor for mask (default 0.9).
IMPORTANT: Notice that the symmetry parameters refer to the surface
itself, and not to the grid.
2. Dependencies of the thorn
This thorn additionally requires thorns: GenericAHF, util.
For testing purposes, one also needs: analyticBH, exact.
3. Thorn distribution
This thorn is available to everyone.
4. Additional information
