add max_allowable_horizon_radius

parameter to (we think) fix a serious bug where if one horizon goes
off the edge of the grid, we fail to find any other horizons


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

5 files changed, 56 insertions, 26 deletions

diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 059265c..51fb8f8 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -70,7 +70,7 @@ void broadcast_horizon_data(const cGH* GH,
 void print_status(int N_active_procs,
 		  const struct iteration_status_buffers& isb);
 void Newton_step(patch_system& ps,
-		 fp max_allowable_Delta_h_over_h,
+		 fp mean_horizon_radius, fp max_allowable_Delta_h_over_h,
 		 const struct verbose_info& verbose_info);
 	  }
 
@@ -183,14 +183,14 @@ if (hs.has_genuine_horizons())
 			      iteration);
 
 	jtutil::norm<fp> Theta_norms;
-	const enum expansion_status
-	  expansion_status = expansion(ps_ptr, add_to_expansion,
-				       cgi, gi,
-				       error_info, (iteration == 1),
-				       true,	// yes, we want Jacobian coeffs
-				       verbose_info.print_algorithm_details,
-				       &Theta_norms);
-	const bool Theta_is_ok = (expansion_status == expansion_success);
+	const enum expansion_status raw_expansion_status
+		= expansion(ps_ptr, add_to_expansion,
+			   cgi, gi,
+			   error_info, (iteration == 1),
+			   true,	// yes, we want Jacobian coeffs
+			   verbose_info.print_algorithm_details,
+			   &Theta_norms);
+	const bool Theta_is_ok = (raw_expansion_status == expansion_success);
 	const bool norms_are_ok = horizon_is_genuine && Theta_is_ok;
 	if (verbose_info.print_algorithm_debug)
 	   then {
@@ -211,7 +211,6 @@ if (hs.has_genuine_horizons())
 			      iteration);
 
 
-
 	//
 	// have we found this horizon?
 	// if so, compute and output BH diagnostics
@@ -241,6 +240,26 @@ if (hs.has_genuine_horizons())
 
 
 	//
+	// compute the mean horizon radius and check if it's too large
+	// ... if so, then pretend expansion() returned a "surface too large"
+	//     error status
+	//
+
+	jtutil::norm<fp> h_norms;
+	if (horizon_is_genuine)
+	   then ps_ptr->ghosted_gridfn_norms(gfns::gfn__h, h_norms);
+	const fp mean_horizon_radius = horizon_is_genuine ? h_norms.mean()
+							  : 0.0;
+	const bool horizon_is_too_large
+		= (mean_horizon_radius > solver_info
+					 .max_allowable_horizon_radius);
+
+	const enum expansion_status effective_expansion_status
+		= horizon_is_too_large ? expansion_failure__surface_too_large
+				       : raw_expansion_status;
+
+
+	//
 	// see if we need more iterations (either on this or another horizon)
 	//
 
@@ -249,6 +268,7 @@ if (hs.has_genuine_horizons())
 	const bool this_horizon_needs_more_iterations
 	   = horizon_is_genuine && Theta_is_ok
 	     && !found_this_horizon
+	     && !horizon_is_too_large
 	     && (iteration < max_iterations);
 
 	// do I (this processor) need to do more iterations
@@ -280,7 +300,7 @@ if (hs.has_genuine_horizons())
 	  = broadcast_status(GH,
 			     N_procs, N_active_procs,
 			     my_proc, my_active_flag,
-			     hn, iteration, expansion_status,
+			     hn, iteration, effective_expansion_status,
 			     (norms_are_ok ? Theta_norms.rms_norm() : 0.0),
 			     (norms_are_ok ? Theta_norms.infinity_norm() : 0.0),
 			     found_this_horizon, I_need_more_iterations,
@@ -467,7 +487,8 @@ if (hs.has_genuine_horizons())
 	// take the Newton step (scaled if need be)
 	//
 	Newton_step(*ps_ptr,
-		    solver_info.max_allowable_Delta_h_over_h,
+		    mean_horizon_radius, solver_info
+					 .max_allowable_Delta_h_over_h,
 		    verbose_info);
 
 	// end of this Newton iteration
@@ -499,7 +520,7 @@ assert( false );
 // N_procs = The total number of processors.
 // N_active_procs = The number of active processors.
 // my_active_flag = Is this processor an active processor?
-// hn,iteration,status,
+// hn,iteration,effective_expansion_status,
 // {rms,infinity}_norm,
 // found_this_horizon,I_need_more_iterations
 //	= On an active processors, these are the values to be broadcast.
@@ -517,7 +538,7 @@ bool broadcast_status(const cGH* GH,
 		      int N_procs, int N_active_procs,
 		      int my_proc, bool my_active_flag,
 		      int hn, int iteration,
-		      enum expansion_status expansion_status,
+		      enum expansion_status effective_expansion_status,
 		      fp rms_norm, fp infinity_norm,
 		      bool found_this_horizon, bool I_need_more_iterations,
 		      struct iteration_status_buffers& isb)
@@ -544,9 +565,9 @@ assert( my_proc < N_procs );
 //
 
 // the buffers are actually 2-D arrays; these are the column numbers
-// ... if we wanted to, we could pack hn, iteration, and expansion_status
-//     all into a single (64-bit) floating-point value, but it's not
-//     worth the trouble...
+// ... if we wanted to, we could pack hn, iteration, and
+//     effective_expansion_status all into a single (64-bit)
+//     floating-point value, but it's not worth the trouble...
 enum	{
 	// CCTK_INT buffer
 	buffer_var__hn = 0,	// also encodes found_this_horizon flag
@@ -596,7 +617,7 @@ if (my_active_flag)
 	send_buffer(my_proc, buffer_var__iteration)
 		= I_need_more_iterations ? +iteration : -iteration;
 	send_buffer(my_proc, buffer_var__expansion_status)
-		= int(expansion_status);
+		= int(effective_expansion_status);
 	send_buffer(my_proc, buffer_var__rms_norm)      = rms_norm;
 	send_buffer(my_proc, buffer_var__infinity_norm) = infinity_norm;
 	}
@@ -851,6 +872,7 @@ void print_status(int N_active_procs,
 //
 // Arguments:
 // ps = The patch system containing the gridfns h and Delta_h.
+// mean_horizon_radius = ||h||_mean
 // max_allowable_Delta_h_over_h = The maximum allowable
 //				     ||Delta_h||_infinity / ||h||_mean
 //				  Any step over this is internally clamped
@@ -858,16 +880,15 @@ void print_status(int N_active_procs,
 //
 namespace {
 void Newton_step(patch_system& ps,
-		 fp max_allowable_Delta_h_over_h,
+		 fp mean_horizon_radius, fp max_allowable_Delta_h_over_h,
 		 const struct verbose_info& verbose_info)
 {
 //
 // compute scale factor (1 for small steps, <1 for large steps)
 //
 
-jtutil::norm<fp> h_norms;
-ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
-const fp max_allowable_Delta_h = max_allowable_Delta_h_over_h * h_norms.mean();
+const fp max_allowable_Delta_h
+	= max_allowable_Delta_h_over_h * mean_horizon_radius;
 
 jtutil::norm<fp> Delta_h_norms;
 ps.gridfn_norms(gfns::gfn__Delta_h, Delta_h_norms);
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 146b49f..c29e733 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -113,6 +113,7 @@ struct	solver_info
 	int max_Newton_iterations__initial,
 	    max_Newton_iterations__subsequent;
 	fp max_allowable_Delta_h_over_h;
+	fp max_allowable_horizon_radius;
 	fp Theta_norm_for_convergence;
 	};
 
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 4058bbe..88c0159 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -288,6 +288,7 @@ solver_info.max_Newton_iterations__initial
 solver_info.max_Newton_iterations__subsequent
 	= max_Newton_iterations__subsequent;
 solver_info.max_allowable_Delta_h_over_h = max_allowable_Delta_h_over_h;
+solver_info.max_allowable_horizon_radius = max_allowable_horizon_radius;
 solver_info.Theta_norm_for_convergence   = Theta_norm_for_convergence;
 
 
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 11ab58e..599d704 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -32,10 +32,7 @@ enum	Jacobian_compute_method
 	};
 
 //
-// this enum describes the status of an
-//	expansion()
-// and/or
-//	expansion_Jacobian()
+// this enum describes the status of an  expansion()  or  expansion_Jacobian()
 // computation which returns (i.e. which doesn't abort the Cactus run)
 //
 enum	expansion_status
@@ -47,6 +44,13 @@ enum	expansion_status
 	// non-finite() (i.e. +/-infinity or NaN) values found in h
 	expansion_failure__surface_nonfinite,
 
+	// surface is too large
+	// ... this value is never returned by  expansion()  or
+	//      expansion_Jacobian() , but is placed in this enum
+	//     for the convenience of higher-level software which
+	//     wishes to have such an error condition
+	expansion_failure__surface_too_large,
+
 	// geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE
 	expansion_failure__surface_outside_grid,
 
diff --git a/src/gr/misc-gr.cc b/src/gr/misc-gr.cc
index b17a7e9..9ca3d1f 100644
--- a/src/gr/misc-gr.cc
+++ b/src/gr/misc-gr.cc
@@ -127,6 +127,9 @@ case expansion_success:
 case expansion_failure__surface_nonfinite:
 	return "infinity/NaN in surface shape!";
 	break;
+case expansion_failure__surface_too_large:
+	return "surface too large";
+	break;
 case expansion_failure__surface_outside_grid:
 	return "surface outside grid";
 	break;