* major changes to driver routines to find multiple horizons in parallel

  across multiple processors -- see src/driver/README.parallel for details
* drop convergence checks on ||Delta_h|| in param.ccl because they don't
  fit well with parallelization changes

==> With this changes, AHFinderDirect is now (I think)
    multiprocessor-ready!!


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

1 files changed, 198 insertions, 464 deletions

diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 01e28e2..fd35afd 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -11,11 +11,7 @@
 /// find_horizon - find a horizon
 /// do_evaluate_expansion
 /// do_test_expansion_Jacobian
-/// do_find_horizon
 ///
-/// compute_BH_diagnostics - compute BH diagnostics for a horizon
-/// surface_integral - compute surface integral of a gridfn over the horizon
-//
 
 #include <stdio.h>
 #include <assert.h>
@@ -48,6 +44,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -68,28 +65,25 @@ const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
 					const char gridfn_name[],
 					bool check_for_NULL = true);
 
-bool do_evaluate_expansion
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps,
-	 bool active_flag, int hn, int N_horizons);
-bool do_test_expansion_Jacobian
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info,
-	 bool test_all_Jacobian_compute_methods,
-	 struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons);
-bool do_find_horizon
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 const struct solver_info& solver_info, bool initial_find_flag,
-	 const struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons);
+void do_evaluate_expansions(int my_proc, int N_horizons,
+			    horizon_sequence& hs,
+			       struct AH_info* const my_AH_info[],
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct IO_info& IO_info,
+			    bool output_h, bool output_Theta,
+			    const struct verbose_info& verbose_info,
+			    int timer_handle);
+void do_test_expansion_Jacobians(int my_proc, int N_horizons,
+				 struct AH_info* const my_AH_info[],
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				       struct Jacobian_info& Jac_info,
+				 bool test_all_Jacobian_compute_methods,
+				 const struct IO_info& IO_info,
+				 const struct verbose_info& verbose_info,
+				 int timer_handle);
+
 
 void compute_BH_diagnostics
 	(const patch_system& ps,
@@ -112,27 +106,30 @@ extern "C"
 DECLARE_CCTK_ARGUMENTS
 DECLARE_CCTK_PARAMETERS
 
-const struct verbose_info& verbose_info = state.verbose_info;
-      struct IO_info&      IO_info      = state.IO_info;
-const struct solver_info&  solver_info  = state.solver_info;
-
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
 
+const int my_proc = state.my_proc;
+horizon_sequence& hs = *state.my_hs;
+const bool active_flag = hs.has_genuine_horizons();
+
+      struct cactus_grid_info&          cgi =  state.cgi;
+const struct    geometry_info&           gi =  state.gi;
+      struct    Jacobian_info&     Jac_info =  state.Jac_info;
+      struct          IO_info&      IO_info =  state.IO_info;
+const struct     verbose_info& verbose_info =  state.verbose_info;
+
 // what are the semantics of the Cactus gxx variables? (these may
 // change from one call to another, so we have to re-check each time)
 if      (CCTK_Equals(metric_type, "physical"))
-   then state.cgi.Cactus_conformal_metric = false;
+   then cgi.Cactus_conformal_metric = false;
 else if (CCTK_Equals(metric_type, "static conformal"))
-   then state.cgi.Cactus_conformal_metric = (conformal_state > 0);
+   then cgi.Cactus_conformal_metric = (conformal_state > 0);
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
 		   metric_type);				/*NOTREACHED*/
 
-// only processor #0 is "active"
-const bool active_flag = (state.my_proc == 0);
-
-// get the Cactus time step and decide if we want to output [hH] now
+// get the Cactus time step and decide if we want to output h and/or Theta now
 IO_info.time_iteration = cctk_iteration;
 IO_info.time           = cctk_time;
 const bool output_h
@@ -142,32 +139,22 @@ const bool output_Theta
    = (IO_info.how_often_to_output_Theta > 0)
      && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0);
 
-//
 // if we're using them,
 //    we need to re-fetch the Cactus data pointers at least each time step,
 //    because they change each time Cactus rotates the time levels
-//
-if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
-   then setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
+if (gi.geometry_method == geometry__local_interp_from_Cactus_grid)
+   then setup_Cactus_gridfn_data_ptrs(cctkGH, cgi);
 
-	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
+// set initial guess for any (genuine) horizons that need it,
+// i.e. for any (genuine) horizons where we didn't find the horizon previously
+	for (int hn = hs.init_hn() ; hs.is_genuine() ; hn = hs.next_hn())
 	{
-	struct AH_info& AH_info = state.AH_info_array[hn];
-	patch_system& ps = *AH_info.ps_ptr;
-
-	//
-	// If this is our first attempt to find this horizon, or
-	//    if we've tried to find it before but we failed on our
-	//    immediately previous attempt, then we need to (re)set
-	//    the initial guess.
-	// Otherwise (i.e. if we've tried to find this horizon before,
-	//    and we succeeded on our immediately previous attempt),
-	//    then we can just reuse the previous horizon position as
-	//    our initial guess for this time around.
-	//
-	const bool initial_find_flag = ! AH_info.AH_found;
-	if (initial_find_flag)
-	   then {
+	assert( state.my_AH_info[hn] != NULL );
+	struct AH_info& AH_info = *state.my_AH_info[hn];
+	if (AH_info.found_flag)
+	   then AH_info.initial_find_flag = false;
+	   else {
+		patch_system& ps = *AH_info.ps_ptr;
 		setup_initial_guess(ps,
 				    AH_info.initial_guess_info,
 				    IO_info,
@@ -177,129 +164,51 @@ if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
 				      IO_info, IO_info.h_base_file_name,
 				      hn, verbose_info
 					  .print_algorithm_highlights);
+		AH_info.initial_find_flag = true;
 		}
+	}
 
-	//
-	// create the Jacobian matrix if we're going to need it
-	// and it doesn't already exist
-	//
-	switch	(state.method)
-		{
-	case method__evaluate_expansion:
-		// no Jacobian needed ==> no-op here
-		break;
-	case method__test_expansion_Jacobian:
-	case method__find_horizon:
-		if (AH_info.Jac_ptr == NULL)
-		   then AH_info.Jac_ptr
-			   = new_Jacobian(state.Jac_info
-					       .Jacobian_store_solve_method,
-					  ps,
-					  verbose_info
-					  .print_algorithm_highlights);
-		break;
-	default:
-		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "\n"
-			   "   find_horizons(): unknown method=(int)%d!\n"
-			   "                    in create-Jacobian switch\n"
-			   "                    (this should never happen!)"
-			   ,
-			   int(state.method));			/*NOTREACHED*/
-		}
+//
+// now the main horizon finding (or other computation)
+//
+switch	(state.method)
+	{
+case method__evaluate_expansions:
+	do_evaluate_expansions(my_proc, state.N_horizons,
+			       *state.my_hs, state.my_AH_info,
+			       cgi, gi,
+			       IO_info, output_h, output_Theta,
+			       state.verbose_info, state.timer_handle);
+	break;
 
-	AH_info.AH_found = false;
-	switch	(state.method)
-		{
-	case method__evaluate_expansion:
-		do_evaluate_expansion(verbose_info, state.timer_handle,
-				      IO_info, output_h, output_Theta,
-				      state.cgi, state.gi,
-				      ps,
-				      active_flag, hn, state.N_horizons);
-		break;
-
-	case method__test_expansion_Jacobian:
-		do_test_expansion_Jacobian(verbose_info, state.timer_handle,
-					   IO_info,
-				       (test_all_Jacobian_compute_methods != 0),
-					   state.Jac_info, state.cgi, state.gi,
-					   ps, *AH_info.Jac_ptr,
-					   active_flag, hn, state.N_horizons);
-		break;
-
-	case method__find_horizon:
-		AH_info.AH_found
-		   = do_find_horizon(verbose_info, state.timer_handle,
-				     IO_info, output_h, output_Theta,
-				     solver_info, initial_find_flag,
-				     state.Jac_info, state.cgi, state.gi,
-				     ps, *AH_info.Jac_ptr,
-				     active_flag, hn, state.N_horizons);
-
-		// store found flag in Cactus array variable
-		AH_found[hn] = AH_info.AH_found;
-
-		if (AH_info.AH_found)
-		   then {
-// compute BH diagnostics
-compute_BH_diagnostics(ps,
-		       state.BH_diagnostics_info,
-		       verbose_info, AH_info.BH_diagnostics);
-const struct BH_diagnostics& BH_diagnostics = AH_info.BH_diagnostics;
-
-			// print BH diagnostics?
-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));
-	}
+case method__test_expansion_Jacobians:
+	do_test_expansion_Jacobians(my_proc, state.N_horizons,
+				    state.my_AH_info,
+				    cgi, gi, Jac_info, 
+				    (test_all_Jacobian_compute_methods != 0),
+				    IO_info, verbose_info, state.timer_handle);
+	break;
 
-// store BH diagnostics in Cactus array variables
-centroid_x[hn]    = BH_diagnostics.centroid_x;
-centroid_y[hn]    = BH_diagnostics.centroid_y;
-centroid_z[hn]    = BH_diagnostics.centroid_z;
-area[hn]          = BH_diagnostics.area;
-m_irreducible[hn] = BH_diagnostics.m_irreducible;
+case method__find_horizons:
+	if (state.timer_handle >= 0)
+	   then CCTK_TimerStartI(state.timer_handle);
+	Newton(cctkGH,
+	       state.N_procs, state.N_active_procs, my_proc,
+	       *state.my_hs, state.my_AH_info,
+	       cgi, gi, Jac_info, state.solver_info,
+	       IO_info, state.BH_diagnostics_info,
+	       verbose_info);
+	if (state.timer_handle >= 0)
+	   then CCTK_TimerStopI(state.timer_handle);
+	break;
 
-// output BH diagnostics?
-if (active_flag && IO_info.output_BH_diagnostics)
-   then {
-	if (AH_info.BH_diagnostics_fileptr == NULL)
-	   then AH_info.BH_diagnostics_fileptr
-		   = setup_BH_diagnostics_output_file(IO_info,
-						      hn, state.N_horizons);
-	output_BH_diagnostics_fn(BH_diagnostics,
-				 IO_info, hn,
-				 AH_info.BH_diagnostics_fileptr);
-	}
-			}
-		   else {
-			if (verbose_info.print_physics_details)
-			   then CCTK_VInfo(CCTK_THORNSTRING,
-					   "no apparent horizon found");
-			}
-		break;
-
-	default:
-		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "\n"
-			   "   find_horizons(): unknown method=(int)%d!\n"
-			   "                    in main switch\n"
-			   "                    (this should never happen!)"
-			   ,
-			   int(state.method));			/*NOTREACHED*/
-		}
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   find_horizons(): unknown method=(int)%d!\n"
+"                    (this should never happen!)"
+		   ,
+		   int(state.method));				/*NOTREACHED*/
 	}
 
 if (state.timer_handle >= 0)
@@ -387,8 +296,10 @@ return data_ptr;
 //******************************************************************************
 
 //
-// This function implements  AHFinderDirect::method == "horizon function",
-// that is, it evaluates the Theta(h) function for a single apparent horizon.
+// This function implements  AHFinderDirect::method == "horizon function":
+// On processor #0 it evaluates the Theta(h) function for each apparent
+// horizon (and does any I/O desired); on other processors it does N_horizons
+// dummy evaluations on horizon #0.
 //
 // Note that if we decide to output h, we output it *after* any Theta(h)
 // evaluation or horizon finding has been done, to ensure that all the
@@ -398,63 +309,85 @@ return data_ptr;
 // timer_handle = a valid Cactus timer handle if we want to time the
 //		  apparent horizon process, or -ve to skip this
 //		  (we only time the computation, not the file I/O)
-// output_[hH] = flags to control whether or not we should output [hH]
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We evaluate Theta(h) as usual, but don't write
-//			   any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the evaluation was successful, false
-// otherwise.
 //
 namespace {
-bool do_evaluate_expansion
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps,
-	 bool active_flag, int hn, int N_horizons)
+void do_evaluate_expansions(int my_proc, int N_horizons,
+			    horizon_sequence& hs,
+			       struct AH_info* const my_AH_info[],
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct IO_info& IO_info,
+			    bool output_h, bool output_Theta,
+			    const struct verbose_info& verbose_info,
+			    int timer_handle)
 {
-jtutil::norm<fp> Theta_norms;
-
-if (timer_handle >= 0)
-   then CCTK_TimerStartI(timer_handle);
-const bool status = expansion(ps, cgi, gi, false, true, &Theta_norms);
-if (timer_handle >= 0)
-   then CCTK_TimerStopI(timer_handle);
-if (!status)
-   then return false;					// *** ERROR RETURN ***
-
-if (Theta_norms.is_nonempty())	// might be empty if Theta(h) eval failed
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
-		   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
-
-if (active_flag && output_h)
-   then output_gridfn(ps, gfns::gfn__h,
-		      IO_info, IO_info.h_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-if (active_flag && output_Theta)
-   then output_gridfn(ps, gfns::gfn__Theta,
-		      IO_info, IO_info.Theta_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-
-return true;						// *** NORMAL RETURN ***
+const bool active_flag = (my_proc == 0);
+
+if (active_flag)
+   then {
+	assert( hs.N_horizons() == N_horizons );
+	assert( hs.my_N_horizons() == N_horizons );
+
+		for (int hn = hs.init_hn() ;
+		     hs.is_genuine() ;
+		     hn = hs.next_hn())
+		{
+		assert( my_AH_info[hn] != NULL );
+		struct AH_info& AH_info = *my_AH_info[hn];
+		patch_system& ps = *AH_info.ps_ptr;
+
+		if (timer_handle >= 0)
+		   then CCTK_TimerStartI(timer_handle);
+		jtutil::norm<fp> Theta_norms;
+		const bool Theta_ok = expansion(&ps,
+						cgi, gi,
+						false,	// no Jacobian coeffs
+						true,	// yes, print msgs
+						&Theta_norms);
+		if (timer_handle >= 0)
+		   then CCTK_TimerStopI(timer_handle);
+
+		if (output_h)
+		   then output_gridfn(ps, gfns::gfn__h,
+				      IO_info, IO_info.h_base_file_name,
+				      hn, verbose_info.print_algorithm_details);
+
+		if (Theta_ok)
+		   then {
+			CCTK_VInfo(CCTK_THORNSTRING,
+			   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
+			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
+			if (output_Theta)
+			   then output_gridfn(ps, gfns::gfn__Theta,
+					      IO_info, IO_info
+						       .Theta_base_file_name,
+					      hn, verbose_info
+						  .print_algorithm_details);
+			}
+		}
+	}
+   else {
+		for (int i = 0 ; i < N_horizons ; ++i)
+		{
+		expansion(NULL,	// dummy computation
+			  cgi, gi);
+		}
+	}
 }
 	  }
 
 //******************************************************************************
 
 //
-// This function implements  AHFinderDirect::method == "test Jacobian",
-// that is, it computes and prints the Jacobian matrix J[Theta(h)] for a
-// single apparent horizon.  The computation may optionally be done
-// in several different ways, in which case all the resulting Jacobian
-// matrices are printed, as are their differences.  Alternatively, only
+// This function implements
+//  AHFinderDirect::method == "test expansion Jacobians":
+// On processor #0 it computes and prints the Jacobian matrix J[Theta(h)]
+// function for horizon #1; on other processors it does dummy Jacobian
+// computations.
+//
+// The Jacobian computation may optionally be done in several different
+// ways, in which case all the resulting Jacobian matrices are printed,
+// as are their differences.  Alternatively, only
 // the numerical perturbation computation may be done/printed.
 //
 // Arguments:
@@ -468,259 +401,60 @@ return true;						// *** NORMAL RETURN ***
 //	  false ==> Just test/print the numerical perturbation calculation.
 //		    (This may be useful if one or more of the other methods
 //		    is broken.) 
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We compute the Jacobian(s) as usual, but don't
-//			   write any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the Jacobian computation was successful,
-// false otherwise.
 //
 namespace {
-bool do_test_expansion_Jacobian
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info,
-	 bool test_all_Jacobian_compute_methods,
-	 struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons)
+void do_test_expansion_Jacobians(int my_proc, int N_horizons,
+				 struct AH_info* const my_AH_info[],
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				       struct Jacobian_info& Jac_info,
+				 bool test_all_Jacobian_compute_methods,
+				 const struct IO_info& IO_info,
+				 const struct verbose_info& verbose_info,
+				 int timer_handle)
 {
-// numerical perturbation
-Jacobian* Jac_NP_ptr = & Jac;
-if (!  expansion(ps, cgi, gi, true))
-   then return false;					// *** ERROR RETURN ***
+const bool active_flag = (my_proc == 0);
+assert(N_horizons >= 1);
+
+const bool print_msg_flag = true;
+const int hn = 1;
+
+struct AH_info* AH_info_ptr = active_flag ? my_AH_info[hn]      : NULL;
+patch_system*        ps_ptr = active_flag ? AH_info_ptr->ps_ptr : NULL;
+
+//
+// numerical-perturbation Jacobian
+//
+Jacobian* Jac_NP_ptr = active_flag ? AH_info_ptr->Jac_ptr : NULL;
+expansion(ps_ptr,
+	  cgi, gi);
 Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
-if (! expansion_Jacobian(ps, *Jac_NP_ptr,
-			 cgi, gi, Jac_info,
-			 true))	// print msgs
-   then return false;					// *** ERROR RETURN ***
+expansion_Jacobian(ps_ptr, Jac_NP_ptr,
+		   cgi, gi, Jac_info,
+		   print_msg_flag);
 
 Jacobian* Jac_SD_FDdr_ptr = NULL;
 if (test_all_Jacobian_compute_methods)
    then {
 	// symbolic differentiation with finite diff d/dr
-	Jac_SD_FDdr_ptr = new_Jacobian(Jac_info.Jacobian_store_solve_method,
-				       ps,
-				       verbose_info.print_algorithm_details);
-	if (! expansion(ps, cgi, gi, true))
-	   then return false;				// *** ERROR RETURN ***
+	Jac_SD_FDdr_ptr = active_flag
+			  ? new_Jacobian(Jac_info.Jacobian_store_solve_method,
+					 *ps_ptr,
+					 verbose_info.print_algorithm_details)
+			  : NULL;
+	expansion(ps_ptr,
+		  cgi, gi,
+		  true);		// compute SD Jacobian coeffs
 	Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
-	if (! expansion_Jacobian(ps, *Jac_SD_FDdr_ptr,
-				 cgi, gi, Jac_info,
-				 true))	// print msgs
-	   then return false;				// *** ERROR RETURN ***
+	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
+			   cgi, gi, Jac_info,
+			   print_msg_flag);
 	}
 
 if (active_flag)
-   then output_Jacobians(ps,
+   then output_Jacobians(*ps_ptr,
 			 Jac_NP_ptr, Jac_SD_FDdr_ptr,
 			 IO_info, IO_info.Jacobian_base_file_name,
-			 hn, true);
-
-return true;						// *** NORMAL RETURN ***
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function implements  AHFinderDirect::method == "find horizon",
-// that is, it finds the (an) apparent horizon.
-//
-// Arguments:
-// timer_handle = a valid Cactus timer handle if we want to time the
-//		  apparent horizon process, or -ve to skip this
-//		  (we only time the computation, not the file I/O)
-// initial_find_flag = true ==> This is the first time we've tried to find
-//				this horizon, or we've tried before but
-//				failed on our most recent previous attempt.
-//				Thus we should use "initial" parameters
-//				for the Newton iteration.
-//		       flase ==> We've tried to find this horizon before,
-//				 and we succeeded on our most recent
-//				 previous attempt.  Thus we should use
-//				 "subsequent" parameters for the Newton
-//				 iteration.
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We find the horizon as usual, but don't write
-//			   any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the apparent horizon was found
-// successfully, false otherwise.
-//
-namespace {
-bool do_find_horizon
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 const struct solver_info& solver_info, bool initial_find_flag,
-	 const struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons)
-{
-if (verbose_info.print_algorithm_highlights)
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "searching for horizon %d/%d",
-		   hn, N_horizons);
-
-if (timer_handle >= 0)
-   then CCTK_TimerStartI(timer_handle);
-const bool status = Newton_solve(ps, Jac,
-				 cgi, gi, Jac_info,
-				 solver_info, initial_find_flag, IO_info,
-				 active_flag, hn, verbose_info);
-if (timer_handle >= 0)
-   then CCTK_TimerStopI(timer_handle);
-if (! status)
-   then return false;					// *** ERROR RETURN ***
-
-if (active_flag && output_h)
-   then output_gridfn(ps, gfns::gfn__h,
-		      IO_info, IO_info.h_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-if (active_flag && output_Theta)
-   then output_gridfn(ps, gfns::gfn__Theta,
-		      IO_info, IO_info.Theta_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-
-return true;						// *** NORMAL RETURN ***
-}
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// Given that an apparent horizon has been found, this function computes
-// various black hole diagnostics.
-//
-// Inputs (gridfns)
-// h		# ghosted
-// one		# nominal
-// global_[xyz]	# nominal
-//
-namespace {
-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 raw 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);
-
-//
-// adjust integrals to take into account patch system not covering full sphere
-//
-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(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n"
-"                             (this should never happen!)"
-,
-		   int(ps.type()));				/*NOTREACHED*/
-	}
-
-BH_diagnostics.centroid_x = integral_x / integral_one;
-BH_diagnostics.centroid_y = integral_y / integral_one;
-BH_diagnostics.centroid_z = integral_z / integral_one;
-
-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));
-}
-	  }
-
-//******************************************************************************
-
-//
-// 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)
-{
-return ps.integrate_gridfn
-	   (src_gfn,
-	    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);
+			 hn, print_msg_flag);
 }
 	  }