compute min/max horizon radius for each horizon found

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

13 files changed, 96 insertions, 102 deletions

diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc
index 5ec1418..9754fe9 100644
--- a/src/driver/BH_diagnostics.cc
+++ b/src/driver/BH_diagnostics.cc
@@ -9,7 +9,7 @@
 // print_BH_diagnostics - print a summary of BH diagnostics
 //
 // setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file
-// output_BH_diagnostics_fn - append BH diagnostics to file
+// do_output_BH_diagnostics - append BH diagnostics to file
 //
 
 #include <stdio.h>
@@ -44,6 +44,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -66,7 +67,7 @@ fp surface_integral(const patch_system& ps, int src_gfn,
 //
 BH_diagnostics::BH_diagnostics()
 	: centroid_x(0.0), centroid_y(0.0), centroid_z(0.0),
-	  mean_radius(0.0),
+	  min_radius(0.0), max_radius(0.0), mean_radius(0.0),
 	  circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0),
 	  area(0.0),
 	  m_irreducible(0.0)
@@ -88,10 +89,17 @@ BH_diagnostics::BH_diagnostics()
 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)
 {
 //
+// compute min/max range 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();
+
+//
 // compute raw surface integrals
 //
 fp integral_one = surface_integral(ps, gfns::gfn__one,
@@ -206,24 +214,20 @@ return ps.integrate_gridfn
 // This function prints a summary of the BH diagnostics, using CCTK_VInfo().
 //
 void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			  int N_horizons, int hn,
-			  const struct verbose_info& verbose_info)
+			  int N_horizons, int hn)
 {
-if (verbose_info.print_physics_details)
-   then {
-	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));
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
-		   hn, N_horizons,
-		   double(BH_diagnostics.area),
-		   double(BH_diagnostics.m_irreducible));
-	}
+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));
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "AH %d/%d: area=%.10g m_irreducible=%.10g",
+	   hn, N_horizons,
+	   double(BH_diagnostics.area),
+	   double(BH_diagnostics.m_irreducible));
 }
 
 //******************************************************************************
@@ -257,17 +261,19 @@ if (fileptr == NULL)
 
 fprintf(fileptr, "# apparent horizon %d/%d\n", hn, N_horizons);
 fprintf(fileptr, "#\n");
-fprintf(fileptr, "# column 1 = cctk_iteration\n");
-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 = mean radius\n");
-fprintf(fileptr, "# column 7 = xy-plane circumference\n");
-fprintf(fileptr, "# column 8 = xz-plane circumference\n");
-fprintf(fileptr, "# column 9 = yz-plane circumference\n");
-fprintf(fileptr, "# column 10 = area\n");
-fprintf(fileptr, "# column 11 = m_irreducible\n");
+fprintf(fileptr, "# column  1 = cctk_iteration\n");
+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 = min radius\n");
+fprintf(fileptr, "# column  7 = max radius\n");
+fprintf(fileptr, "# column  8 = mean radius\n");
+fprintf(fileptr, "# column  9 = xy-plane circumference\n");
+fprintf(fileptr, "# column 10 = xz-plane circumference\n");
+fprintf(fileptr, "# column 11 = yz-plane circumference\n");
+fprintf(fileptr, "# column 12 = area\n");
+fprintf(fileptr, "# column 13 = m_irreducible\n");
 fflush(fileptr);
 
 return fileptr;
@@ -288,21 +294,22 @@ return fileptr;
 // BH_diagnostics = The BH diagnostics to be written
 // fileptr = The stdio file pointer to append to
 //
-void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
+void do_output_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
 			      const struct IO_info& IO_info,
 			      int hn, FILE* fileptr)
 {
 assert(fileptr != NULL);
 
 fprintf(fileptr,
-//  cctk_iteration        mean radius                     area    m_irreducible
-//  ==  cctk_time         ======  {xy,xz,yz}-plane        ======  ======
-//  ==  ====  centroid_[xyz]      circumferences          ======  ======
-//  ==  ====  ==========  ======  ======================  ======  ======
-   "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
+//  cctk_iteration        mean radius     max radius                      area    m_irreducible
+//  ==  cctk_time         ======  min radius      {xy,xz,yz}-plane        ======  ======
+//  ==  ====  centroid_[xyz]      ======  ======  circumferences          ======  ======
+//  ==  ====  ==========  ======  ======  ======  ======================  ======  ======
+   "%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,
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 27ec3b0..5283a44 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -39,6 +39,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -287,13 +288,14 @@ if (hs.has_genuine_horizons())
 	if (found_this_horizon)
 	   then {
 		compute_BH_diagnostics(*ps_ptr, BH_diagnostics_info,
-				       verbose_info,
 				       AH_info_ptr->BH_diagnostics);
-		// 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,
-				     verbose_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);
+			}
 		if (IO_info.output_BH_diagnostics)
 		   then {
 			if (AH_info_ptr->BH_diagnostics_fileptr == NULL)
@@ -301,7 +303,7 @@ if (hs.has_genuine_horizons())
 				  = setup_BH_diagnostics_output_file(IO_info,
 								     N_horizons,
 								     hn);
-			output_BH_diagnostics_fn(AH_info_ptr->BH_diagnostics,
+			do_output_BH_diagnostics(AH_info_ptr->BH_diagnostics,
 						 IO_info,
 						 hn, AH_info_ptr
 						     ->BH_diagnostics_fileptr);
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 225130a..0995f50 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -4,6 +4,8 @@
 //
 // prerequisites:
 //	<stdio.h>
+//	"horizon_sequence.hh"
+//	"BH_diagnostics.hh"
 //
 
 //******************************************************************************
@@ -93,14 +95,6 @@ struct	solver_info
 	fp Theta_norm_for_convergence;
 	};
 
-//
-// This struct holds info for computing black hole diagnostics.
-//
-struct	BH_diagnostics_info
-	{
-	enum patch::integration_method integral_method;
-	};
-
 //******************************************************************************
 
 //
@@ -164,23 +158,6 @@ struct	verbose_info
 //******************************************************************************
 
 //
-// (A single copy of) this struct 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
-	{
-	fp centroid_x, centroid_y, centroid_z;
-	fp mean_radius;
-	fp circumference_xy, circumference_xz, circumference_yz;
-	fp area;
-	fp m_irreducible;
-
-	// constructor initializes all diagnostics to 0.0
-	BH_diagnostics();
-	};
-
-//
 // (A single copy of) this struct holds all of our information about
 // a single apparent horizon.
 //
@@ -284,21 +261,6 @@ void Newton(cGH* GH,
 	    const struct error_info& error_info,
 	    const struct verbose_info& verbose_info);
 
-// BH_diagnostics.cc
-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);
-void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
-			  int N_horizons, int hn,
-			  const struct verbose_info& verbose_info);
-FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
-				       int N_horizons, int hn);
-void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
-			      const struct IO_info& IO_info,
-			      int hn, FILE* fileptr);
-
 // io.cc
 void input_gridfn(patch_system& ps, int unknown_gfn,
 		  const struct IO_info& IO_info, const char base_file_name[],
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 0324cf1..5db1698 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -45,6 +45,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
index 7e00390..9f563dd 100644
--- a/src/driver/initial_guess.cc
+++ b/src/driver/initial_guess.cc
@@ -42,6 +42,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/driver/io.cc b/src/driver/io.cc
index 3db7edf..2a68725 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -5,7 +5,7 @@
 // output_gridfn - write an angular grid function to an output file
 // output_Jacobians - write a Jacobian matrix or matrices to an output file
 //
-/// io_file_name
+/// io_file_name - compute file name for angular-gridfn I/O file
 //
 
 #include <stdio.h>
@@ -41,6 +41,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -301,7 +302,7 @@ fclose(fileptr);
 //******************************************************************************
 
 //
-// This function encapsulates our file-naming conventions for "horizon"
+// This function encapsulates our file-naming conventions for angular-gridfn
 // output files (those used for h, H, and other angular grid functions).
 //
 // Arguments:
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index d0b823e..77cdd10 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -48,6 +48,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/driver/state.cc b/src/driver/state.cc
index a6092c7..b7e64ae 100644
--- a/src/driver/state.cc
+++ b/src/driver/state.cc
@@ -33,6 +33,7 @@ using jtutil::error_exit;
 #include "../gr/gr.hh"
 
 #include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/jtutil/norm.cc b/src/jtutil/norm.cc
index edd3fd1..6fd3a92 100644
--- a/src/jtutil/norm.cc
+++ b/src/jtutil/norm.cc
@@ -19,10 +19,15 @@ namespace jtutil
 template <typename fp_t>
   void norm<fp_t>::data(fp_t x)
 {
-++N_;
 sum_  += x;
 sum2_ += x*x;
-infinity_norm_ = jtutil::max(infinity_norm_, jtutil::abs<fp_t>(x));
+
+const fp_t abs_x = jtutil::abs<fp_t>(x);
+max_abs_value_ = jtutil::max(max_abs_value_, abs_x);
+min_abs_value_ = (N_ == 0) ? abs_x
+			   : jtutil::min(min_abs_value_, abs_x);
+
+++N_;
 }
 	  }	// namespace jtutil::
 
@@ -38,8 +43,6 @@ template<typename fp_t>
 template<typename fp_t>
   fp_t norm<fp_t>::rms_norm() const
 	{ assert(is_nonempty()); return sqrt(sum2_/fp_t(N_)); }
-template<typename fp_t>
-  fp_t norm<fp_t>::infinity_norm() const { return infinity_norm_; }
 	  }	// namespace jtutil::
 
 //******************************************************************************
diff --git a/src/jtutil/test_norm.cc b/src/jtutil/test_norm.cc
index b156873..25a9c59 100644
--- a/src/jtutil/test_norm.cc
+++ b/src/jtutil/test_norm.cc
@@ -25,6 +25,8 @@ foo.data(-3.0);
 assert( fuzzy<double>::EQ(foo.two_norm(), sqrt(45.0)) );
 assert( fuzzy<double>::EQ(foo.rms_norm(), 3.0) );
 assert( fuzzy<double>::EQ(foo.infinity_norm(), 3.0) );
+assert( fuzzy<double>::EQ(foo.max_abs_value(), 3.0) );
+assert( fuzzy<double>::EQ(foo.min_abs_value(), 3.0) );
 
 foo.reset();
 assert(  foo.is_empty() );
@@ -40,6 +42,8 @@ foo.data(-5.0);
 assert( fuzzy<double>::EQ(foo.two_norm(), sqrt(52.0)) );
 assert( fuzzy<double>::EQ(foo.rms_norm(), sqrt(10.4)) );
 assert( fuzzy<double>::EQ(foo.infinity_norm(), 5.0) );
+assert( fuzzy<double>::EQ(foo.max_abs_value(), 5.0) );
+assert( fuzzy<double>::EQ(foo.min_abs_value(), 1.0) );
 
 printf("everything ok!\n");
 return 0;
diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh
index d87330d..2012a28 100644
--- a/src/jtutil/util.hh
+++ b/src/jtutil/util.hh
@@ -88,7 +88,8 @@ double modulo_reduce(double x, double xmod, double xmin, double xmax);
 //******************************************************************************
 
 //
-// This template class computes means, 2-norms, rms-norms, and infinity-norms.
+// This template class computes means, 2-norms, rms-norms, and
+// infinity-norms.  It also computes minimum values.
 //
 template <typename fp_t>
 class	norm
@@ -98,7 +99,9 @@ public:
 	fp_t mean() const;
 	fp_t two_norm() const;		// sqrt(sum x_i^2)
 	fp_t rms_norm() const;		// sqrt(average of x_i^2)
-	fp_t infinity_norm() const;	// max(|x_i|)
+	fp_t infinity_norm() const { return max_abs_value_; }
+	fp_t max_abs_value() const { return max_abs_value_; }
+	fp_t min_abs_value() const { return min_abs_value_; }
 
 	// specify data point
 	void data(fp_t x);
@@ -108,12 +111,20 @@ public:
 	bool is_nonempty() const { return N_ > 0; }
 
 	// reset ==> just like newly-constructed object
-	void reset() { N_ = 0; sum2_ = 0.0; infinity_norm_ = 0.0; }
+	void reset()
+		{
+		N_ = 0;
+		sum2_ = 0.0;
+		max_abs_value_ = 0.0;
+		min_abs_value_ = 0.0;
+		}
 
 	// constructor, destructor
 	// ... compiler-generated no-op destructor is ok
 	norm()
-		: N_(0), sum_(0.0), sum2_(0.0), infinity_norm_(0.0)
+		: N_(0),
+		  sum_(0.0), sum2_(0.0),
+		  max_abs_value_(0.0), min_abs_value_(0.0)
 		{ }
 
 private:
@@ -127,7 +138,8 @@ private:
 	int N_;			// # of data points
 	fp_t sum_;		// sum(data)
 	fp_t sum2_;		// sum(data^2)
-	fp_t infinity_norm_;	// max |data|
+	fp_t max_abs_value_;	// max |data|
+	fp_t min_abs_value_;	// min |data|
 	};
 
 //******************************************************************************
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 7046b30..3adaaed 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -1512,7 +1512,8 @@ norms.reset();
 //******************************************************************************
 
 //
-// This function computes norms of a ghosted-grid gridfn.
+// This function computes norms of a ghosted-grid gridfn over the
+// nominal grid.
 //
 void patch_system::ghosted_gridfn_norms(int ghosted_src_gfn,
 					jtutil::norm<fp>& norms)
@@ -1527,12 +1528,10 @@ norms.reset();
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
 	const patch& p = ith_patch(pn);
-		for (int irho = p.ghosted_min_irho() ;
-		     irho <= p.ghosted_max_irho() ;
-		     ++irho)
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
 		{
-		for (int isigma = p.ghosted_min_isigma() ;
-		     isigma <= p.ghosted_max_isigma() ;
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
 		norms.data(p.ghosted_gridfn(ghosted_src_gfn, irho,isigma));
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 99bead8..ce7b100 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -272,7 +272,7 @@ public:
 	// dst += delta
 	void add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn);
 
-	// compute norms of gridfn
+	// compute norms of gridfn (only over nominal grid)
 	void         gridfn_norms(int         src_gfn, jtutil::norm<fp>& norms)
 		const;
 	void ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms)