This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@3 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 397 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..a61721d
--- /dev/null
+++ b/README
@@ -0,0 +1,397 @@
+Cactus Code Thorn AHFinder
+Authors    : Miguel Alcubierre
+Managed by : Miguel Alcubierre <miguel@aei-potsdam.mpg.de>
+Version    : 1.0 
+--------------------------------------------------------------------------
+
+
+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.
+
+
+1. Purpose of the thorn
+
+MAIN 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.
+
+
+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 "mahf_flow" equal to "yes":
+
+mahf_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 lines to the
+parameter file:
+
+ah_persists       = "yes"
+minimumAHF_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:
+
+mahfgrid  =  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 -2.0).
+
+
+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
+