add compuatation and printing of mean radius for each AH found

output format for BH info in stdout log is now

INFO (AHFinderDirect): AH 1/1: r=1.86525 at (-0.000000,0.000000,-0.000000)
INFO (AHFinderDirect): AH 1/1: area=45.23876324 m_irreducible=0.9486815054

and in BH_diagnostics file is now

# apparent horizon 1 of 1
#
# column 1 = cctk_iteration
# column 2 = cctk_time
# column 3 = centroid_x
# column 4 = centroid_y
# column 5 = centroid_z
# column 6 = mean radius
# column 7 = area
# column 8 = m_irreducible


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

8 files changed, 58 insertions, 29 deletions

diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 3db180d..53526a5 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -169,6 +169,7 @@ struct	BH_diagnostics
 	{
 	fp centroid_x, centroid_y, centroid_z;
 	fp area;
+	fp mean_radius;
 	fp m_irreducible;
 	};
 
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index bfeedef..6ff104b 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -219,18 +219,21 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 				= AH_info.BH_diagnostics;
 
 			// print BH diagnostics?
-			if (verbose_info.print_physics_details)
-			   then {
-				CCTK_VInfo(CCTK_THORNSTRING,
-					   "AH found: A=%.10g at (%f,%f,%f)",
-					   double(BH_diagnostics.area),
-					   double(BH_diagnostics.centroid_x),
-					   double(BH_diagnostics.centroid_y),
-					   double(BH_diagnostics.centroid_z));
-				CCTK_VInfo(CCTK_THORNSTRING,
-					   "m_irreducible=%.10g",
-					  double(BH_diagnostics.m_irreducible));
-				}
+if (verbose_info.print_physics_details)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: r=%g at (%f,%f,%f)",
+		   hn, state.N_horizons,
+		   double(BH_diagnostics.mean_radius),
+		   double(BH_diagnostics.centroid_x),
+		   double(BH_diagnostics.centroid_y),
+		   double(BH_diagnostics.centroid_z));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
+		   hn, state.N_horizons,
+		   double(BH_diagnostics.area),
+		   double(BH_diagnostics.m_irreducible));
+	}
 
 			// store BH diagnostics in Cactus array variables
 			centroid_x[hn]    = BH_diagnostics.centroid_x;
@@ -519,7 +522,7 @@ Jacobian& Jac = *Jac_ptr;
 
 if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "searching for horizon %d of %d",
+		   "searching for horizon %d/%d",
 		   hn, N_horizons);
 
 if (timer_handle >= 0)
@@ -572,6 +575,8 @@ void compute_BH_diagnostics
 //
 fp integral_one = surface_integral(ps, gfns::gfn__one,
 				   BH_diagnostics_info.surface_integral_method);
+fp integral_h = surface_integral(ps, gfns::gfn__h,
+				 BH_diagnostics_info.surface_integral_method);
 fp integral_x = surface_integral(ps, gfns::gfn__global_x,
 				 BH_diagnostics_info.surface_integral_method);
 fp integral_y = surface_integral(ps, gfns::gfn__global_y,
@@ -588,6 +593,7 @@ case patch_system::patch_system__full_sphere:
 	break;
 case patch_system::patch_system__plus_z_hemisphere:
 	integral_one *= 2.0;
+	integral_h *= 2.0;
 	integral_x *= 2.0;
 	integral_y *= 2.0;
 	integral_z = integral_one * ps.origin_z();
@@ -595,12 +601,14 @@ case patch_system::patch_system__plus_z_hemisphere:
 case patch_system::patch_system__plus_xy_quadrant_mirrored:
 case patch_system::patch_system__plus_xy_quadrant_rotating:
 	integral_one *= 4.0;
+	integral_h *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
 	integral_z *= 4.0;
 	break;
 case patch_system::patch_system__plus_xz_quadrant_rotating:
 	integral_one *= 4.0;
+	integral_h *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y *= 4.0;
 	integral_z = integral_one * ps.origin_z();
@@ -608,6 +616,7 @@ case patch_system::patch_system__plus_xz_quadrant_rotating:
 case patch_system::patch_system__plus_xyz_octant_mirrored:
 case patch_system::patch_system__plus_xyz_octant_rotating:
 	integral_one *= 8.0;
+	integral_h *= 8.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
 	integral_z = integral_one * ps.origin_z();
@@ -626,6 +635,7 @@ BH_diagnostics.centroid_y = integral_y / integral_one;
 BH_diagnostics.centroid_z = integral_z / integral_one;
 
 BH_diagnostics.area = integral_one;
+BH_diagnostics.mean_radius = integral_h / integral_one;
 BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
 }
 	  }
diff --git a/src/driver/io.cc b/src/driver/io.cc
index 5b4aaaa..1cad171 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -352,8 +352,9 @@ fprintf(fileptr, "# column 2 = cctk_time\n");
 fprintf(fileptr, "# column 3 = centroid_x\n");
 fprintf(fileptr, "# column 4 = centroid_y\n");
 fprintf(fileptr, "# column 5 = centroid_z\n");
-fprintf(fileptr, "# column 6 = area\n");
-fprintf(fileptr, "# column 7 = m_irreducible\n");
+fprintf(fileptr, "# column 6 = mean radius\n");
+fprintf(fileptr, "# column 7 = area\n");
+fprintf(fileptr, "# column 8 = m_irreducible\n");
 fflush(fileptr);
 
 return fileptr;
@@ -383,14 +384,16 @@ assert(fileptr != NULL);
 //               cctk_iteration
 //               ==   cctk_time
 //               ==   ====  centroid_[xyz]
-//               ==   ====  ==========  area   m_irreducible
-//               ==   ====  ==========  =====  =====
-fprintf(fileptr, "%d\t%.3f\t%f\t%f\t%f\t%.10g\t%.10g\n",
+//               ==   ====  ==========  mean_radius
+//               ==   ====  ==========  =====  area   m_irreducible
+//               ==   ====  ==========  =====  =====  ====
+fprintf(fileptr, "%d\t%.3f\t%f\t%f\t%f\t%.10g\t%.10g\t%.10g\n",
 	IO_info.time_iteration, double(IO_info.time),
 	BH_diagnostics.centroid_x,
 	BH_diagnostics.centroid_y,
 	BH_diagnostics.centroid_z,
-	BH_diagnostics.area, BH_diagnostics.m_irreducible);
+	BH_diagnostics.mean_radius, BH_diagnostics.area,
+	BH_diagnostics.m_irreducible);
 
 fflush(fileptr);
 }
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 482e8bb..f15319b 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -253,7 +253,7 @@ state.AH_info_ptrs.push_back(NULL);
 	{
 	if (verbose_info.print_algorithm_highlights)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "   setting up data structures for horizon %d of %d",
+			   "   setting up data structures for horizon %d/%d",
 			   hn, int(N_horizons));
 
 	struct AH_info* AH_ptr = new AH_info;
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index 5067ccb..931adf6 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -269,7 +269,10 @@ else	error_exit(ERROR_EXIT,
 // where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma).
 //
 // Arguments:
-// src_gfn = (in) The gridfn to be integrated.
+// unknown_src_gfn = (in) The gridfn to be integrated.  This may be
+//			  either nominal-grid or ghosted-grid; n.b. in
+//			  the latter case the integral is still done only
+//			  over the patch's nominal area.
 // ghosted_radius_gfn = (in) The surface radius.
 // g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
 // method = (in) Selects the integration scheme.
@@ -278,7 +281,7 @@ else	error_exit(ERROR_EXIT,
 // - The integration coefficients are computed in a rather inefficient
 //   manner.
 //
-fp patch::integrate_gridfn(int src_gfn,
+fp patch::integrate_gridfn(int unknown_src_gfn,
 			   int ghosted_radius_gfn,
 			   int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 					 int g_yy_gfn, int g_yz_gfn,
@@ -286,12 +289,15 @@ fp patch::integrate_gridfn(int src_gfn,
 			   enum integration_method method)
 	const
 {
+const bool src_is_ghosted = is_valid_ghosted_gfn(unknown_src_gfn);
+
 fp sum = 0.0;
 	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
 	{
 	for (int isigma = min_isigma() ; isigma <= max_isigma() ; ++isigma)
 	{
-	const fp fn = gridfn(src_gfn, irho,isigma);
+	const fp fn = unknown_gridfn(src_is_ghosted,
+				     unknown_src_gfn, irho,isigma);
 
 	const fp rho   = rho_of_irho    (irho);
 	const fp sigma = sigma_of_isigma(isigma);
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 621a29f..07958eb 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -362,12 +362,16 @@ public:
 	  enum integration_method
 	    decode_integration_method(const char method_string[]);
 
-	// compute the surface integral of a gridfn over the patch,
+	// compute the surface integral of a gridfn over the patch's
+	// nominal area,
 	//	$\int f(\rho,\sigma) \, dA$
 	//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
 	// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma)
 	// ... integration method selected by  method  argument
-	fp integrate_gridfn(int src_gfn,
+	// ... src gridfn may be either nominal-grid or ghosted-grid
+	//     (n.b. in the latter case the integral is still done
+	//           only over the patch's nominal area)
+	fp integrate_gridfn(int unknown_src_gfn,
 			    int ghosted_radius_gfn,
 			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 					  int g_yy_gfn, int g_yz_gfn,
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index da4027c..0aefef6 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -1560,12 +1560,13 @@ norms.reset();
 // overlap, i.e. that patch_overlap_width=1
 //
 // Arguments:
-// src_gfn = (in) The gridfn to be integrated.
+// unknown_src_gfn = (in) The gridfn to be integrated.  This may be
+//			  either nominal-grid or ghosted-grid.
 // ghosted_radius_gfn = (in) The surface radius.
 // g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
 // method = (in) Selects the integration scheme.
 //
-fp patch_system::integrate_gridfn(int src_gfn,
+fp patch_system::integrate_gridfn(int unknown_src_gfn,
 				  int ghosted_radius_gfn,
 				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 						int g_yy_gfn, int g_yz_gfn,
@@ -1577,7 +1578,7 @@ fp integral = 0.0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
 	const patch& p = ith_patch(pn);
-	integral += p.integrate_gridfn(src_gfn,
+	integral += p.integrate_gridfn(unknown_src_gfn,
 				       ghosted_radius_gfn,
 				       g_xx_gfn, g_xy_gfn, g_xz_gfn,
 						 g_yy_gfn, g_yz_gfn,
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 6bd5518..d4fe3ce 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -283,11 +283,15 @@ public:
 	//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
 	// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma)
 	// ... integration method selected by  method  argument
+	// ... src gridfn may be either nominal-grid or ghosted-grid
 	// ... note integral is only over patch system itself,
 	//     even if this covers a proper subset of the full sphere
 	//     FIXME: it would be better to compute the full-sphere integral
 	//            (but this would require knowing the gridfn's symmetries)
-	fp integrate_gridfn(int src_gfn,
+	//     FIXME: current implementation assumes patches' nominal grids
+	//            form a partition of the region covered, i.e. they
+	//            just touch each other
+	fp integrate_gridfn(int unknown_src_gfn,
 			    int ghosted_radius_gfn,
 			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 					  int g_yy_gfn, int g_yz_gfn,