* make BH_diagnostics more of an object

  (eg move non-member code into member fns)
* handle expansion() and expansion_Jacobian() now returning full status
  instead of just Boolean success/failure flag


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

4 files changed, 171 insertions, 199 deletions

diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc
index 9754fe9..06f8423 100644
--- a/src/driver/BH_diagnostics.cc
+++ b/src/driver/BH_diagnostics.cc
@@ -3,13 +3,12 @@
 //
 // BH_diagnostics::BH_diagnostics - initialize a  struct BH_diagnostics
 //
-// compute_BH_diagnostics - compute BH diagnostics
-/// surface_integral - compute surface integral over a patch system
+// BH_diagnostics::compute - compute BH diagnostics after an AH has been found
+/// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere
 //
-// print_BH_diagnostics - print a summary of BH diagnostics
-//
-// setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file
-// do_output_BH_diagnostics - append BH diagnostics to file
+// print - print a line or two summarizing the diagnostics
+// setup_output_file - create/open output file, write header describing fields
+// output - write a (long) line of all the diagnostics
 //
 
 #include <stdio.h>
@@ -48,17 +47,6 @@ using jtutil::error_exit;
 #include "driver.hh"
 
 //******************************************************************************
-
-//
-// prototypes for functions local to this file
-//
-
-namespace {
-fp surface_integral(const patch_system& ps, int src_gfn,
-		    enum patch::integration_method method);
-	  }
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -86,104 +74,72 @@ BH_diagnostics::BH_diagnostics()
 // one		# nominal
 // global_[xyz]	# nominal
 //
-void compute_BH_diagnostics
+void BH_diagnostics::compute
 	(const patch_system& ps,
-	 const struct BH_diagnostics_info& BH_diagnostics_info,
-	 struct BH_diagnostics& BH_diagnostics)
+	 const struct BH_diagnostics_info& BH_diagnostics_info)
 {
 //
-// compute min/max range of horizon
+// min/max radius of horizon
 //
 jtutil::norm<fp> h_norms;
 ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
-BH_diagnostics.min_radius = h_norms.min_abs_value();
-BH_diagnostics.max_radius = h_norms.max_abs_value();
+min_radius = h_norms.min_abs_value();
+max_radius = h_norms.max_abs_value();
 
 //
-// compute raw surface integrals
+// surface integrals
 //
-fp integral_one = surface_integral(ps, gfns::gfn__one,
-				   BH_diagnostics_info.integral_method);
-fp integral_h = surface_integral(ps, gfns::gfn__h,
-				 BH_diagnostics_info.integral_method);
-fp integral_x = surface_integral(ps, gfns::gfn__global_x,
-				 BH_diagnostics_info.integral_method);
-fp integral_y = surface_integral(ps, gfns::gfn__global_y,
-				 BH_diagnostics_info.integral_method);
-fp integral_z = surface_integral(ps, gfns::gfn__global_z,
-				 BH_diagnostics_info.integral_method);
+const fp integral_one = surface_integral(ps,
+					 gfns::gfn__one, true, true, true,
+					 BH_diagnostics_info.integral_method);
+const fp integral_h = surface_integral(ps,
+				       gfns::gfn__h, true, true, true,
+				       BH_diagnostics_info.integral_method);
+const fp integral_x = surface_integral(ps,
+				       gfns::gfn__global_x, true, true, false,
+				       BH_diagnostics_info.integral_method);
+const fp integral_y = surface_integral(ps,
+				       gfns::gfn__global_y, true, false, true,
+				       BH_diagnostics_info.integral_method);
+const fp integral_z = surface_integral(ps,
+				       gfns::gfn__global_z, false, true, true,
+				       BH_diagnostics_info.integral_method);
 
 //
-// adjust integrals to take into account patch system not covering full sphere
+// centroids
 //
-switch	(ps.type())
-	{
-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();
-	break;
-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();
-	break;
-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();
-	break;
-default:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n"
-"                             (this should never happen!)"
-,
-		   int(ps.type()));				/*NOTREACHED*/
-	}
+centroid_x = integral_x / integral_one;
+centroid_y = integral_y / integral_one;
+centroid_z = integral_z / integral_one;
 
-BH_diagnostics.centroid_x = integral_x / integral_one;
-BH_diagnostics.centroid_y = integral_y / integral_one;
-BH_diagnostics.centroid_z = integral_z / integral_one;
+//
+// area, mean radius, and mass
+//
+area = integral_one;
+mean_radius = integral_h / integral_one;
+m_irreducible = sqrt(area / (16.0*PI));
 
-BH_diagnostics.area = integral_one;
-BH_diagnostics.circumference_xy = ps.xy_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.circumference_xz = ps.xz_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.circumference_yz = ps.yz_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.mean_radius = integral_h / integral_one;
-BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
+//
+// circumferences
+//
+circumference_xy
+  = ps.circumference("xy", gfns::gfn__h,
+		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+							     gfns::gfn__g_dd_33,
+		     BH_diagnostics_info.integral_method);
+circumference_xz
+  = ps.circumference("xz", gfns::gfn__h,
+		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+							     gfns::gfn__g_dd_33,
+		     BH_diagnostics_info.integral_method);
+circumference_yz
+  = ps.circumference("yz", gfns::gfn__h,
+		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+							     gfns::gfn__g_dd_33,
+		     BH_diagnostics_info.integral_method);
 }
 
 //******************************************************************************
@@ -192,42 +148,45 @@ BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
 // This function computes the surface integral of a gridfn over the
 // horizon.
 //
-namespace {
-fp surface_integral(const patch_system& ps, int src_gfn,
-		    enum patch::integration_method method)
+//static
+  fp BH_diagnostics::surface_integral
+	(const patch_system& ps,
+	 int src_gfn, bool src_gfn_is_even_across_xy_plane,
+		      bool src_gfn_is_even_across_xz_plane,
+		      bool src_gfn_is_even_across_yz_plane,
+	 enum patch::integration_method method)
 {
 return ps.integrate_gridfn
-	   (src_gfn,
+	   (src_gfn, src_gfn_is_even_across_xy_plane,
+		     src_gfn_is_even_across_xz_plane,
+		     src_gfn_is_even_across_yz_plane,
 	    gfns::gfn__h,
 	    gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
 				gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
 						    gfns::gfn__g_dd_33,
 	    method);
 }
-	  }
 
 //******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
 //
-// This function prints a summary of the BH diagnostics, using CCTK_VInfo().
+// This function prints a line or two summarizing the diagnostics,
+// using CCTK_VInfo().
 //
-void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			  int N_horizons, int hn)
+void BH_diagnostics::print(int N_horizons, int hn)
+	const
 {
 CCTK_VInfo(CCTK_THORNSTRING,
 	   "AH %d/%d: r=%g at (%f,%f,%f)",
 	   hn, N_horizons,
-	   double(BH_diagnostics.mean_radius),
-	   double(BH_diagnostics.centroid_x),
-	   double(BH_diagnostics.centroid_y),
-	   double(BH_diagnostics.centroid_z));
+	   double(mean_radius),
+	   double(centroid_x), double(centroid_y), double(centroid_z));
 CCTK_VInfo(CCTK_THORNSTRING,
 	   "AH %d/%d: area=%.10g m_irreducible=%.10g",
 	   hn, N_horizons,
-	   double(BH_diagnostics.area),
-	   double(BH_diagnostics.m_irreducible));
+	   double(area), double(m_irreducible));
 }
 
 //******************************************************************************
@@ -236,11 +195,14 @@ CCTK_VInfo(CCTK_THORNSTRING,
 
 //
 // This function creates a BH-diagnostics output file, writes a suitable
-// header comment, and returns a stdio file pointer which can be used to
-// append data to the file.
+// header comment identifying the fields to be written by  output() ,
+// flushes the stream (to help in examining the output while Cactus is
+// still running), and finally returns a stdio file pointer which can be
+// used by  output()  to output data to the file.
 //
-FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
-				       int N_horizons, int hn)
+FILE* BH_diagnostics::setup_output_file(const struct IO_info& IO_info,
+					int N_horizons, int hn)
+	const
 {
 const int N_file_name_buffer = 200;
 char file_name_buffer[N_file_name_buffer];
@@ -253,7 +215,7 @@ FILE *fileptr = fopen(file_name_buffer, "w");
 if (fileptr == NULL)
    then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   setup_BH_diagnostics_output_file():\n"
+"   BH_diagnostics::setup_output_file():\n"
 "        can't open BH-diagnostics output file\n"
 "        \"%s\"!"
 		   ,
@@ -282,21 +244,16 @@ return fileptr;
 //******************************************************************************
 
 //
-// This function appends a BH-diagnostics line to a data file.  It
-// attempts to ensure that the newly-written line is flushed to disk,
-// so the output file can be examined while the Cactus run is still in
-// progress.
-//
-// (The "_fn" in the function name is to distinguish it from the Cactus
-// parameter  output_BH_diagnostics .)
+// This function outputs a BH-diagnostics line to a stdio stream, then
+// flushes the stream (to help in examining the output while Cactus is
+// still running).
 //
 // Arguments:
 // BH_diagnostics = The BH diagnostics to be written
 // fileptr = The stdio file pointer to append to
 //
-void do_output_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			      const struct IO_info& IO_info,
-			      int hn, FILE* fileptr)
+void BH_diagnostics::output(FILE*fileptr, const struct IO_info& IO_info)
+	const
 {
 assert(fileptr != NULL);
 
@@ -307,14 +264,10 @@ fprintf(fileptr,
 //  ==  ====  ==========  ======  ======  ======  ======================  ======  ======
    "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\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.min_radius, BH_diagnostics.max_radius,
-	BH_diagnostics.mean_radius,
-	BH_diagnostics.circumference_xy, BH_diagnostics.circumference_xz,
-					 BH_diagnostics.circumference_yz,
-	BH_diagnostics.area,
-	BH_diagnostics.m_irreducible);
+	centroid_x, centroid_y, centroid_z,
+	min_radius, max_radius, mean_radius,
+	circumference_xy, circumference_xz, circumference_yz,
+	area, m_irreducible);
 
 fflush(fileptr);
 }
diff --git a/src/driver/BH_diagnostics.hh b/src/driver/BH_diagnostics.hh
index 259736a..40bd6c7 100644
--- a/src/driver/BH_diagnostics.hh
+++ b/src/driver/BH_diagnostics.hh
@@ -1,6 +1,11 @@
 // BH_diagnostics.hh -- header file for BH diagnostics
 // $Header$
 
+//
+// prerequisites:
+//	<stdio.h>
+//
+
 //******************************************************************************
 
 //
@@ -14,37 +19,60 @@ struct	BH_diagnostics_info
 //******************************************************************************
 
 //
-// (A single copy of) this struct holds all of our black hole diagnostics
+// A  struct BH_diagnostics  holds all of our black hole diagnostics
 // for a single apparent horizon.  These diagnostics are only meaningful
 // if the apparent horizon has indeed been found.
 //
 struct	BH_diagnostics
 	{
+public:
 	fp centroid_x, centroid_y, centroid_z;
 	fp min_radius, max_radius, mean_radius;
+
 	fp circumference_xy, circumference_xz, circumference_yz;
 	fp area;
 	fp m_irreducible;
 
+public:
+	// compute diagnostics (assuming that apparent horizon has been found)
+	void compute(const patch_system& ps,
+		     const struct BH_diagnostics_info& BH_diagnostics_info);
+
+	// print (CCTK_VInfo()) a line or two summarizing diagnostics
+	void print(int N_horizons, int hn)
+		const;
+
+	// create/open output file and write header describing output() fields
+	// ... stream is flushed after output to help with
+	//     looking at diagnostics while Cactus is still running
+	FILE* setup_output_file(const struct IO_info& IO_info,
+				int N_horizons, int hn)
+		const;
+
+	// output a (long) line of all the diagnostics, to a stdio stream
+	// ... stream is flushed after output to help with
+	//     looking at diagnostics while Cactus is still running
+	void output(FILE* fileptr, const struct IO_info& IO_info)
+		const;
+
 	// constructor initializes all diagnostics to 0.0
 	BH_diagnostics();
-	};
 
-//******************************************************************************
+	// no destructor needed, compiler-generated no-op is fine
 
-//
-// prototypes for non-class functions visible outside their source files
-//
+private:
+	// helper function: compute surface integral of specified gridfn
+	static
+	  fp surface_integral(const patch_system& ps,
+			      int src_gfn, bool src_gfn_is_even_across_xy_plane,
+					   bool src_gfn_is_even_across_xz_plane,
+					   bool src_gfn_is_even_across_yz_plane,
+			      enum patch::integration_method method);
 
-// BH_diagnostics.cc
-void compute_BH_diagnostics
-	(const patch_system& ps,
-	 const struct BH_diagnostics_info& BH_diagnostics_info,
-	 struct BH_diagnostics& BH_diagnostics);
-void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			  int N_horizons, int hn);
-FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
-				       int N_horizons, int hn);
-void do_output_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			      const struct IO_info& IO_info,
-			      int hn, FILE* fileptr);
+private:
+	// we forbid copying and passing by value
+	// by declaring the copy constructor and assignment operator
+	// private, but never defining them
+	BH_diagnostics(const BH_diagnostics& rhs);
+	BH_diagnostics& operator=(const BH_diagnostics& rhs);
+	};
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 5283a44..6486ec6 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -175,6 +175,8 @@ if (hs.has_genuine_horizons())
 	//
 	for (int iteration = 1 ; ; ++iteration)
 	{
+	enum expansion_status status;
+
 	if (verbose_info.print_algorithm_debug)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "beginning iteration %d (horizon_is_genuine=%d)",
@@ -193,13 +195,13 @@ if (hs.has_genuine_horizons())
 			      iteration);
 
 	jtutil::norm<fp> Theta_norms;
-	const bool Theta_is_ok
-	   = expansion(ps_ptr,
-		       cgi, gi,
-		       error_info, (iteration == 1),
-		       true,	// yes, we want the Jacobian coeffs
-		       verbose_info.print_algorithm_details,
-		       &Theta_norms);
+	status = expansion(ps_ptr,
+			   cgi, gi,
+			   error_info, (iteration == 1),
+			   true,	// yes, we want the Jacobian coeffs
+			   verbose_info.print_algorithm_details,
+			   &Theta_norms);
+	const bool Theta_is_ok = (status == expansion_success);
 	if (verbose_info.print_algorithm_debug)
 	   then {
 		CCTK_VInfo(CCTK_THORNSTRING,
@@ -287,26 +289,23 @@ if (hs.has_genuine_horizons())
 				       rms_norm_buffer, infinity_norm_buffer);
 	if (found_this_horizon)
 	   then {
-		compute_BH_diagnostics(*ps_ptr, BH_diagnostics_info,
-				       AH_info_ptr->BH_diagnostics);
+		struct BH_diagnostics& BH_diagnostics
+			= AH_info_ptr->BH_diagnostics;
+		BH_diagnostics.compute(*ps_ptr, BH_diagnostics_info);
 		if (verbose_info.print_physics_details)
 		   then {
 			// FIXME: see header comment -- user probably won't see
 			//        this for my_proc != 0 in a multiprocessor run
-			print_BH_diagnostics(AH_info_ptr->BH_diagnostics,
-					     N_horizons, hn);
+			BH_diagnostics.print(N_horizons, hn);
 			}
 		if (IO_info.output_BH_diagnostics)
 		   then {
 			if (AH_info_ptr->BH_diagnostics_fileptr == NULL)
 			   then AH_info_ptr->BH_diagnostics_fileptr
-				  = setup_BH_diagnostics_output_file(IO_info,
-								     N_horizons,
-								     hn);
-			do_output_BH_diagnostics(AH_info_ptr->BH_diagnostics,
-						 IO_info,
-						 hn, AH_info_ptr
-						     ->BH_diagnostics_fileptr);
+				  = BH_diagnostics.setup_output_file
+					(IO_info, N_horizons, hn);
+			BH_diagnostics.output(AH_info_ptr->BH_diagnostics_fileptr,
+					      IO_info);
 			}
 
 		if (IO_info.output_h)
@@ -350,14 +349,14 @@ if (hs.has_genuine_horizons())
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "   computing Jacobian: genuine/dummy flag %d",
 			   int(this_horizon_needs_more_iterations));
-	const bool Jacobian_is_ok
-	   = expansion_Jacobian(this_horizon_needs_more_iterations ? ps_ptr
-								   : NULL,
-				this_horizon_needs_more_iterations ? Jac_ptr
-								   : NULL,
-				cgi, gi, Jacobian_info,
-				error_info, (iteration == 1),
-				verbose_info.print_algorithm_details);
+	status = expansion_Jacobian(this_horizon_needs_more_iterations ? ps_ptr
+								       : NULL,
+				    this_horizon_needs_more_iterations ? Jac_ptr
+								       : NULL,
+				    cgi, gi, Jacobian_info,
+				    error_info, (iteration == 1),
+				    verbose_info.print_algorithm_details);
+	const bool Jacobian_is_ok = (status == expansion_success);
 
 
 	//
@@ -476,18 +475,19 @@ assert( my_proc < N_procs );
 //
 // We do the combined reduce/broadcast in a single reduction operation
 // by setting up a 2-D buffer whose entries are all 0s, except that on each
-// processor the [my_proc] row hase the values we want to reduce/broadcast,
-// then doing a sum-reduction of the buffers across processors.
+// processor the [my_proc] row has the values we want to reduce/broadcast,
+// then doing a sum-reduction of the buffers across processors.  Alas,
+// the input and output buffers must be distinct, so for the broadcast
+// we need two copies of the buffer on each processor.
 //
-// Alas, to do everything in a single operation, all the values have to
-// be of a single datatype in a single array, so we convert  hn  and
-//  iteration  to CCTK_REAL, and encode the Boolean flags in other
-// variables' signs.  Also, for purposes of the reduction we treat the
-// buffer as a 1-D array.
+// To do everything in a single reduction operation, all the values have
+// to be of a single datatype in a single array, so we convert  hn  and
+//  iteration  to CCTK_REAL, and encode the Boolean flags in some of the
+// variables' signs.  Also, for simplicity we treat the buffer as a 1-D
+// array for the reduction.
 //
 
 // the buffer is actually a 2-D array; these are the column numbers
-// ... n.b. input and output buffers must be distinct!
 enum	{
 	buffer_var__hn = 0,
 	buffer_var__iteration,
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 5db1698..869c5d3 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -85,15 +85,6 @@ void do_test_expansion_Jacobians(int my_proc, int N_horizons,
 				 const struct error_info& error_info,
 				 const struct verbose_info& verbose_info,
 				 int timer_handle);
-
-
-void compute_BH_diagnostics
-	(const patch_system& ps,
-	 const struct BH_diagnostics_info& BH_diagnostics_info,
-	 const struct verbose_info& verbose_info,
-	 struct BH_diagnostics& BH_diagnostics);
-fp surface_integral(const patch_system& ps, int src_gfn,
-		    enum patch::integration_method method);
 	  }
 
 //******************************************************************************