add support for computing BH diagnostics (mass, area, centroid)

via surface integrals over the horizon


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@715 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 94 insertions, 30 deletions

diff --git a/param.ccl b/param.ccl
index f899a1c..c7701c0 100644
--- a/param.ccl
+++ b/param.ccl
@@ -3,6 +3,14 @@
 
 ################################################################################
 
+# we may need to look at grid::domain to choose our patch system symmetries
+shares: grid
+USES KEYWORD domain
+
+# all remaining parameters are private to this thorn
+private:
+################################################################################
+
 #
 # overall parameters for the apparent horizon finding algorithm itself
 #
@@ -19,22 +27,27 @@ boolean find_AHs_at_poststep \
 
 keyword method "top-level method used to find the apparent horizon"
 {
-# options mostly for testing/debugging
+# these options are mostly for testing/debugging
 "horizon function"  :: "evaluate the LHS function H(h)"
 "Jacobian test"	    :: \
   "compute/print the J[H(h)] Jacobian matrix by all possible methods"
 "Jacobian test (NP only)" :: \
   "compute/print the J[H(h)] Jacobian matrix by numerical perturbation only"
-# normal AH finding
+
+# this is for normal apparent horizon finding
 "Newton solve"	    :: "find the horizon via Newton's method"
 } "Newton solve"
 
 #
-# at present we support up to 4 horizons
+# At present we support up to 4 horizons; for user-friendliness we
+# number these 1, 2, 3, ... .   Since Cactus arrays are 0-origin, we
+# make the arrays be of size N_horizons+1, and don't use the [0] array
+# elements.
 #
-# to increase this, just raise the upper limit for this variable,
-# and change all occurences of "[4]" in this file to whatever the new
-# upper limit is (no changes are needed in the source code)
+# To change the N_horizons limit, just change the upper limit for
+# N_horizons itself, change all the [N_horizons+1] array sizes in this
+# paramter file, and recompile your configuration.  No changes are
+# needed to the source code.
 #
 int N_horizons "number of apparent horizons to search for"
 {
@@ -51,14 +64,22 @@ int N_horizons "number of apparent horizons to search for"
 keyword verbose_level \
   "controls which (how many) messages to print describing AH finding"
 {
+# 1 line each time step giving number of horizons found and their masses
 "physics highlights"   :: "just a few physics messages"
+
+# 1 line for each horizon giving position/mass/area, + a summary line or two
 "physics details"      :: "more detailed physics messages"
+
+# 1 line giving H(h) norms at each Newton iteration
 "algorithm highlights" :: \
   "physics details + a few messages about the AH-finding algorithm"
+
+# lots of details tracing what the code is doing
 "algorithm details"    :: \
   "physics details + lots of messages about the AH-finding algorithm"
 } "physics details"
 
+# n.b. printing timing stats is independent of  verbose_level
 boolean print_timing_stats \
   "should we print timing stats for the whole apparent-horizon-finding process?"
 {
@@ -76,6 +97,8 @@ keyword Jacobian_method "how do we compute the Jacobian matrix?"
   "*very* slow, but useful for debugging"
 "symbolic differentiation with finite diff d/dr" :: \
   "fast, tricky programming, uses only gij, dx gij, Kij"
+
+# alas, this isn't implemented yet :(
 "symbolic differentiation"  :: \
   "fast, tricky programming, uses gij, dx gij, dxx gij, Kij, dx Kij"
 } "symbolic differentiation with finite diff d/dr"
@@ -191,21 +214,23 @@ string Jacobian_base_file_name "base file name for Jacobian output file(s)"
 #
 private:
 
-real origin_x[4] "global x coordinate of patch system origin"
+real origin_x[5] "global x coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
-real origin_y[4] "global y coordinate of patch system origin"
+real origin_y[5] "global y coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
-real origin_z[4] "global z coordinate of patch system origin"
+real origin_z[5] "global z coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
 
 keyword patch_system_type "what type of patch system should we use?"
 {
+"match Cactus grid symmetry" :: \
+  "choose automagically based on grid symmetries and the patch system's origin"
 "full sphere"	:: "full sphere, no symmetries"
 "+z hemisphere"	:: "mirror symmetry across z=0 plane"
 "+xy quadrant"	:: "90 degree periodic rotation symmetry about z axis"
@@ -213,7 +238,7 @@ keyword patch_system_type "what type of patch system should we use?"
 		    180 degree periodic rotation symmetry about z axis"
 "+xyz octant"	:: "mirror symmetry across z=0 plane *and* \
 		    90 degree periodic rotation symmetry about z axis"
-} "full sphere"
+} "match Cactus grid symmetry"
 
 int N_ghost_points "number of ghost zones on each side of a patch"
 {
@@ -226,10 +251,38 @@ int N_overlap_points \
 1:*:2 :: "any integer >= 0; current implementation requires that it be odd"
 } 1
 
+#
+# In practice the error in the horizon position is usually dominated
+# by the errors from interpolating the Cactus gij and Kij to the horizon
+# position, not by the angular finite differencing or interpatch interpolation
+# errors.  Thus this parameter can be made quite large (low resolution)
+# for better performance, without seriously affecting the accuracy
+# with which we can locate the horizon.
+#
 real delta_drho_dsigma "angular grid spacing of patches, in degrees"
 {
 (0.0:* :: "any real number > 0.0"
-} 5.0
+} 9.0
+
+################################################################################
+
+#
+# parameters for how we compute surface integrals over the horizon
+#
+
+# ... N is the number of grid zones in a patch, in either direction
+keyword surface_integral_method \
+  "how do we compute surface integrals over the horizon?"
+{
+"trapezoid"		:: ""
+"trapezoid rule"	:: "trapezoid rule (2nd order for smooth functions)"
+"Simpson"		:: ""
+"Simpson's rule"	:: \
+  "Simpson's rule (4th order for smooth fns, requires N to be even)"
+"Simpson (variant)"	:: ""
+"Simpson's rule (variant)":: \
+  "Simpson's rule variant (4th order for smooth fns, requires N >= 7)"
+} "trapezoid"
 
 ################################################################################
 
@@ -242,6 +295,7 @@ keyword geometry_method "how do we compute the slice's geometry?"
 {
 # for normal use
 "interpolate from Cactus grid" :: "interpolate gij and Kij from Cactus grid"
+
 # for testing/debugging
 "Schwarzschild/EF"             :: \
   "hard-wire to Schwarzschild spacetime / Eddington-Finkelstein slice"
@@ -252,6 +306,10 @@ keyword geometry_method "how do we compute the slice's geometry?"
 #
 # parameters for geometry_method = "interpolate from Cactus grid"
 #
+# Note that the interpolated gij and Kij should be at least C1, otherwise
+# the Newton iteration may fail to converge all the way down to tight
+# error tolerances.  In practice a Hermite interpolant works well.
+#
 
 string geometry_interpolator_name \
   "name under which the geometry interpolation operator is registered in Cactus"
@@ -320,7 +378,12 @@ real geometry__Schwarzschild_EF__Delta_xyz \
 ################################################################################
 
 #
-# parameters for the interpatch interpolator
+# parameters for the (1-D angular) interpatch interpolator
+#
+# Note there's no necessary relationship between this interpolator and
+# the geometry interpolator.  In particular, this interpolator could
+# reasonably use a different interpolation operator and/or order from
+# the geometry interpolator.
 #
 string interpatch_interpolator_name \
   "name under which the interpatch interpolation operator is registered in Cactus"
@@ -344,10 +407,11 @@ keyword initial_guess_method \
   "method used to set up initial guess for apparent horizon shape"
 {
 "read from file"   :: "read from input file"
-"sphere"           :: "set up a coordinate sphere"
+"Kerr/Kerr"	   :: "set to the horizon of Kerr spacetime in Kerr coordinates"
+"Kerr/Kerr-Schild" :: \
+  "set to the horizon of Kerr spacetime in Kerr-Schild coordinates"
 "ellipsoid"	   :: "set to a coordinate ellipsoid"
-"Kerr/Kerr"	   :: "set to horizon of Kerr spacetime in Kerr coords"
-"Kerr/Kerr-Schild" :: "set to horizon of Kerr spacetime in Kerr-Schild coords"
+"sphere"           :: "set to a coordinate sphere"
 } "read from file"
 
 boolean output_initial_guess \
@@ -356,69 +420,69 @@ boolean output_initial_guess \
 } "true"
 
 # parameters for initial_guess_method = "sphere"
-real initial_guess__sphere__x_center[4] "x coordinate of sphere center"
+real initial_guess__sphere__x_center[5] "x coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__y_center[4] "y coordinate of sphere center"
+real initial_guess__sphere__y_center[5] "y coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__z_center[4] "z coordinate of sphere center"
+real initial_guess__sphere__z_center[5] "z coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__radius[4] "radius of sphere"
+real initial_guess__sphere__radius[5] "radius of sphere"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
 
 # parameters for initial_guess_method = "ellipsoid"
-real initial_guess__ellipsoid__x_center[4] "x coordinate of ellipsoid center"
+real initial_guess__ellipsoid__x_center[5] "x coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__y_center[4] "y coordinate of ellipsoid center"
+real initial_guess__ellipsoid__y_center[5] "y coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__z_center[4] "z coordinate of ellipsoid center"
+real initial_guess__ellipsoid__z_center[5] "z coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__x_radius[4] "x radius of ellipsoid"
+real initial_guess__ellipsoid__x_radius[5] "x radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
-real initial_guess__ellipsoid__y_radius[4] "y radius of ellipsoid"
+real initial_guess__ellipsoid__y_radius[5] "y radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
-real initial_guess__ellipsoid__z_radius[4] "z radius of ellipsoid"
+real initial_guess__ellipsoid__z_radius[5] "z radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
 
 # parameters for initial_guess_method = "Kerr/Kerr" and "Kerr/Kerr-Schild"
-real initial_guess__Kerr__x_posn[4] "x coordinate of Kerr BH"
+real initial_guess__Kerr__x_posn[5] "x coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__Kerr__y_posn[4] "y coordinate of Kerr BH"
+real initial_guess__Kerr__y_posn[5] "y coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__Kerr__z_posn[4] "z coordinate of Kerr BH"
+real initial_guess__Kerr__z_posn[5] "z coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
 # n.b. my convention is that a=J/m^2 is dimensionless,
 #      while MTW take a=J/m=m * (my a)
-real initial_guess__Kerr__mass[4] "mass of Kerr BH"
+real initial_guess__Kerr__mass[5] "mass of Kerr BH"
 {
 (0.0:*		:: "BH mass = any real number > 0"
 } 1.0
-real initial_guess__Kerr__spin[4] "dimensionless spin J/m^2 of Kerr BH"
+real initial_guess__Kerr__spin[5] "dimensionless spin J/m^2 of Kerr BH"
 {
 (-1.0:1.0)	:: "BH spin = J/m^2 = any real number with absolute value < 1"
 } 0.6