expansion() and expansion_Jacobian() now return an enum expansion_status

instead of just a Boolean success/fail flag --> this will let higher-level
code know why a failing evaluation failed, so it can potentially print
a message to tell the user that (eg) the ghost size is too small


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

3 files changed, 152 insertions, 123 deletions

diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index 1946cfa..b64a0b4 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -56,11 +56,12 @@ using jtutil::pow4;
 
 namespace {
 void setup_xyz_posns(patch_system& ps, bool print_msg_flag);
-bool interpolate_geometry(patch_system* ps_ptr,
-			  const struct cactus_grid_info& cgi,
-			  const struct geometry_info& gi,
-			  const struct error_info& error_info, bool initial_flag,
-			  bool print_msg_flag);
+enum expansion_status
+  interpolate_geometry(patch_system* ps_ptr,
+		       const struct cactus_grid_info& cgi,
+		       const struct geometry_info& gi,
+		       const struct error_info& error_info, bool initial_flag,
+		       bool print_msg_flag);
 void convert_conformal_to_physical(patch_system& ps,
 				   bool print_msg_flag);
 void Schwarzschild_EF_geometry(patch_system& ps,
@@ -142,23 +143,17 @@ bool compute_Theta(patch_system& ps, bool Jacobian_flag,
 //			   to compute various (gridwise) norms of Theta.
 //
 // Results:
-// This function returns
-// * true for a successful computation
-// * true for a dummy computation (ps_ptr == NULL)
-// * false if the computation failed for any of the following reasons:
-//   - the geometry interpolation would need data outside the Cactus grid,
-//     or from an excised region
-//     FIXME: excision isn't implemented yet :(
-//   - NaNs are found in h or in the interpolate geometry variables
-//   - Theta_D <= 0 (this means the interpolated g_ij aren't positive definite)
+// This function returns a status code indicating whether the computation
+// succeeded or failed, and if the latter, what caused the failure.
 //
-bool expansion(patch_system* ps_ptr,
-	       const struct cactus_grid_info& cgi,
-	       const struct geometry_info& gi,
-	       const struct error_info& error_info, bool initial_flag,
-	       bool Jacobian_flag /* = false */,
-	       bool print_msg_flag /* = false */,
-	       jtutil::norm<fp>* Theta_norms_ptr /* = NULL */)
+enum expansion_status
+  expansion(patch_system* ps_ptr,
+	    const struct cactus_grid_info& cgi,
+	    const struct geometry_info& gi,
+	    const struct error_info& error_info, bool initial_flag,
+	    bool Jacobian_flag /* = false */,
+	    bool print_msg_flag /* = false */,
+	    jtutil::norm<fp>* Theta_norms_ptr /* = NULL */)
 {
 const bool active_flag = (ps_ptr != NULL);
 if (print_msg_flag)
@@ -178,7 +173,8 @@ if (active_flag)
 	if (gi.check_that_h_is_finite && !h_is_finite(*ps_ptr,
 						      error_info, initial_flag,
 						      print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
+	   then return expansion_failure__surface_nonfinite;
+							// *** ERROR RETURN ***
 
 	// set up xyz positions of grid points
 	setup_xyz_posns(*ps_ptr, print_msg_flag);
@@ -191,11 +187,16 @@ case geometry__global_interp_from_Cactus_grid:
 case geometry__local_interp_from_Cactus_grid:
 	// this is the only function we call unconditionally; it looks at
 	// ps_ptr (non-NULL vs NULL) to choose a normal vs dummy computation
-	if (!interpolate_geometry(ps_ptr,
-				  cgi, gi,
-				  error_info, initial_flag,
-				  print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
+	  {
+	const
+	  enum expansion_status status
+	    = interpolate_geometry(ps_ptr,
+				   cgi, gi,
+				   error_info, initial_flag,
+				   print_msg_flag);
+	if (status != expansion_success)
+	   then return status;				// *** ERROR RETURN ***
+	  }
 	if (active_flag && cgi.Cactus_conformal_metric)
 	   then convert_conformal_to_physical(*ps_ptr,
 					      print_msg_flag);
@@ -222,7 +223,8 @@ if (active_flag)
 	    && !geometry_is_finite(*ps_ptr,
 				   error_info, initial_flag,
 				   print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
+	   then return expansion_failure__geometry_nonfinite;
+							// *** ERROR RETURN ***
 
 	// compute remaining gridfns --> $\Theta$
 	// and optionally also the Jacobian coefficients
@@ -231,10 +233,11 @@ if (active_flag)
 			   Jacobian_flag, Theta_norms_ptr,
 			   error_info, initial_flag,
 			   print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
+	   then return expansion_failure__gij_not_positive_definite;
+							// *** ERROR RETURN ***
 	}
 
-return true;						// *** NORMAL RETURN ***
+return expansion_success;				// *** NORMAL RETURN ***
 }
 
 //******************************************************************************
@@ -328,17 +331,15 @@ if (print_msg_flag)
 // This function may also modify the interpolator parameter table.
 //
 // Results:
-// This function returns
-// * true for a successful computation
-// * true for a dummy computation (ps_ptr == NULL)
-// * false if the computation failed for any of the following reasons:
-//   - the geometry interpolation would need data outside the Cactus grid,
-//     or from an excised region
-//     FIXME: excision isn't implemented yet :(
-// Any other interpolator error return codes are treated as a fatal error.
+// This function returns a status code indicating whether the computation
+// succeeded or failed, and if the latter, what caused the failure.  Possible
+// failure codes are
+// * expansion_failure__surface_outside_grid
+// * expansion_failure__surface_in_excised_region	// not implemented yet
 //
 namespace {
-bool interpolate_geometry(patch_system* ps_ptr,
+enum expansion_status
+  interpolate_geometry(patch_system* ps_ptr,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  const struct error_info& error_info, bool initial_flag,
@@ -692,7 +693,7 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE)
 			   global_x, global_y, global_z,
 			   direction, axis);
 		}
-	return false;					// *** ERROR RETURN ***
+	return expansion_failure__surface_outside_grid;	// *** ERROR RETURN ***
 	}
 if (status < 0)
    then error_exit(ERROR_EXIT,
@@ -701,7 +702,7 @@ if (status < 0)
 ,
 		   status);					/*NOTREACHED*/
 
-return true;						// *** NORMAL RETURN ***
+return expansion_success;				// *** NORMAL RETURN ***
 }
 	  }
 
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index ef2f7e8..7d102f3 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -48,7 +48,8 @@ using jtutil::error_exit;
 //
 
 namespace {
-bool expansion_Jacobian_NP
+enum expansion_status
+  expansion_Jacobian_NP
 	(patch_system& ps, Jacobian& Jac,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
@@ -68,7 +69,8 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
 			     const patch& xp, const ghost_zone& xmgz,
 			     int x_II,
 			     int xm_irho, int xm_isigma);
-bool expansion_Jacobian_dr_FD
+enum expansion_status
+  expansion_Jacobian_dr_FD
 	(patch_system* ps_ptr, Jacobian* Jac_ptr,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
@@ -99,32 +101,32 @@ bool expansion_Jacobian_dr_FD
 // Jac_ptr --> The Jacobian, or == NULL to do (only) a dummy computation.
 //
 // Results:
-// This function returns
-// * true for a successful computation
-// * true for a dummy computation (ps_ptr == NULL)
-// * false if the computation failed for any of the following reasons
-//   documented for  expansion()  (in "expansion.cc").
-//
-bool expansion_Jacobian(patch_system* ps_ptr,
-			Jacobian* Jac_ptr,
-			const struct cactus_grid_info& cgi,
-			const struct geometry_info& gi,
-			const struct Jacobian_info& Jacobian_info,
-			const struct error_info& error_info, bool initial_flag,
-			bool print_msg_flag /* = false */)
+// This function returns a status code indicating whether the computation
+// succeeded or failed, and if the latter, what caused the failure.
+//
+enum expansion_status
+  expansion_Jacobian(patch_system* ps_ptr,
+		     Jacobian* Jac_ptr,
+		     const struct cactus_grid_info& cgi,
+		     const struct geometry_info& gi,
+		     const struct Jacobian_info& Jacobian_info,
+		     const struct error_info& error_info, bool initial_flag,
+		     bool print_msg_flag /* = false */)
 {
 const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL);
+enum expansion_status status;
 
 switch	(Jacobian_info.Jacobian_compute_method)
 	{
 case Jacobian__numerical_perturbation:
 	if (active_flag)
 	   then {
-		if (! expansion_Jacobian_NP(*ps_ptr, *Jac_ptr,
-					    cgi, gi, Jacobian_info,
-					    error_info, initial_flag,
-					    print_msg_flag))
-		   then return false;			// *** ERROR RETURN ***
+		status = expansion_Jacobian_NP(*ps_ptr, *Jac_ptr,
+					       cgi, gi, Jacobian_info,
+					       error_info, initial_flag,
+					       print_msg_flag);
+		if (status != expansion_success)
+		   then return status;			// *** ERROR RETURN ***
 		break;
 		}
 	   else error_exit(ERROR_EXIT,
@@ -137,8 +139,7 @@ case Jacobian__symbolic_diff:
 	error_exit(ERROR_EXIT,
 "***** expansion_Jacobian():\n"
 "        Jacobian_compute_method == \"symbolic differentiation\"\n"
-"        isn't implemented (yet); reverting to\n"
-"        \"symbolic differentiation with finite diff d/dr\"\n");/*NOTREACHED*/
+"        isn't implemented (yet)!\n");				/*NOTREACHED*/
 
 case Jacobian__symbolic_diff_with_FD_dr:
 	if (active_flag)
@@ -147,11 +148,14 @@ case Jacobian__symbolic_diff_with_FD_dr:
 					      print_msg_flag);
 	// this function looks at ps_ptr and Jac_ptr (non-NULL vs NULL)
 	// to choose a normal vs dummy computation
-	if (! expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr,
-				       cgi, gi, Jacobian_info,
-				       error_info, initial_flag,
-				       print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
+	  {
+	status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr,
+					  cgi, gi, Jacobian_info,
+					  error_info, initial_flag,
+					  print_msg_flag);
+	if (status != expansion_success)
+	   then return status;				// *** ERROR RETURN ***
+	  }
 	break;
 
 default:
@@ -163,7 +167,7 @@ default:
 		   int(Jacobian_info.Jacobian_compute_method));	/*NOTREACHED*/
 	}
 
-return true;						// *** NORMAL RETURN ***
+return expansion_success;				// *** NORMAL RETURN ***
 }
 
 //******************************************************************************
@@ -200,12 +204,12 @@ return true;						// *** NORMAL RETURN ***
 //	As implied by the above algorithm, it's traversed by columns.
 //
 // Results:
-// This function returns true for a successful computation, or false
-// for failure.  See  expansion()  (in "expansion.cc")
-// for details of the various ways the computation may fail.
+// This function returns a status code indicating whether the computation
+// succeeded or failed, and if the latter, what caused the failure.
 //
 namespace {
-bool expansion_Jacobian_NP
+enum expansion_status
+  expansion_Jacobian_NP
 	(patch_system& ps, Jacobian& Jac,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
@@ -238,10 +242,12 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 
 	const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
 	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
-	if (! expansion(&ps,
-			cgi, gi,
-			error_info, initial_flag))
-	   then return false;				// *** ERROR RETURN ***
+	const
+	  enum expansion_status status = expansion(&ps,
+						   cgi, gi,
+						   error_info, initial_flag);
+	if (status != expansion_success)
+	   then return status;				// *** ERROR RETURN ***
 
 		for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
 		{
@@ -271,7 +277,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
    	} 
 
 ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
-return true;						// *** NORMAL RETURN ***
+return expansion_success;				// *** NORMAL RETURN ***
 }
 	  }
 
@@ -520,16 +526,12 @@ patch& yp = ye.my_patch();
 //	Jac += d/dr terms
 //
 // Results:
-// This function returns
-// * true for a successful computation
-// * true for a dummy computation (ps_ptr == NULL)
-// * false if the computation failed for any of the following reasons:
-//   - the geometry interpolation would need data outside the Cactus grid,
-//     or from an excised region
-//     FIXME: excision isn't implemented yet :(
+// This function returns a status code indicating whether the computation
+// succeeded or failed, and if the latter, what caused the failure.
 //
 namespace {
-bool expansion_Jacobian_dr_FD
+enum expansion_status
+  expansion_Jacobian_dr_FD
 	(patch_system* ps_ptr, Jacobian* Jac_ptr,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
@@ -551,10 +553,12 @@ if (active_flag)
 	ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 	ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
 	}
-if (! expansion(ps_ptr,
-		cgi, gi,
-		error_info, initial_flag))
-   then return false;					// *** ERROR RETURN ***
+const
+  enum expansion_status status = expansion(ps_ptr,
+					   cgi, gi,
+					   error_info, initial_flag);
+if (status != expansion_success)
+   then return status;					// *** ERROR RETURN ***
 
 if (active_flag)
    then {
@@ -583,6 +587,6 @@ if (active_flag)
 	ps_ptr->gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
 	}
 
-return true;						// *** NORMAL RETURN ***
+return expansion_success;				// *** NORMAL RETURN ***
 }
 	  }
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 58275f5..daf233e 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -31,6 +31,36 @@ enum	Jacobian_compute_method
 	Jacobian__symbolic_diff // no comma
 	};
 
+//
+// this enum describes the status of an
+//	expansion()
+// and/or
+//	expansion_Jacobian()
+// computation which returns (i.e. which doesn't abort the Cactus run)
+//
+enum	expansion_status
+	{
+	// successful computation
+	// ... this is also returned for dummy computations
+	expansion_success,
+
+	// non-finite() (i.e. +/-infinity or NaN) values found in h
+	expansion_failure__surface_nonfinite,
+
+	// geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE
+	expansion_failure__surface_outside_grid,
+
+	// geometry interpolator returns CCTK_ERROR_INTERP_POINT_EXCISED
+	expansion_failure__surface_in_excised_region,	// (not implemented yet)
+
+	// non-finite (i.e. +/-infinity or NaN) values found
+	// in interpolated geometry (g_ij, K_ij, partial_k g_ij)
+	expansion_failure__geometry_nonfinite,
+
+	// interpolated g_ij isn't positive definite
+	expansion_failure__gij_not_positive_definite // no comma
+	};
+
 //******************************************************************************
 
 //
@@ -142,26 +172,26 @@ struct	Jacobian_info
 //
 struct	error_info
 	{
-	// at present, we "handle" all these errors/warnings by CCTK_VWarn();
-	// these parameters give the warning levels for various conditions:
-
-	// point outside (or too close to boundary of) Cactus grid,
-	// i.e. geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE
+	// CCTK_VWarn() level for point outside (or too close to any
+	// boundary of) the Cactus grid, i.e. geometry interpolator returns
+	// CCTK_ERROR_INTERP_POINT_OUTSIDE
 	// ... warning level if error occurs on first Newton iteration,
 	//     i.e. when evaluating expansion/Jacobian for initial guess
 	int warn_level__point_outside__initial;
 	// ... warning level if error occurs on a subsequent Newton iteration
 	int warn_level__point_outside__subsequent;
 
-	// the Cactus configure process didn't find a finite() function,
-	// so we're skipping the "check that the geometry is finite" tests
-	// even though the user set  check_that_geometry_is_finite = "true"
+	// CCTK_VWarn() level for the user set
+	//  check_that_geometry_is_finite = "true"
+	// but the Cactus configure process failed to find the finite()
+	// function ==> we have to skip the check
 	int warn_level__skipping_finite_check;
 
-	// infinity or NaN in interpolated {g_ij, partial_k g_ij, K_ij}
+	// CCTK_VWarn() level for infinity or NaN in interpolated geometry
+	// (g_ij, partial_k g_ij, and/or K_ij)
 	int warn_level__nonfinite_geometry;
 
-	// interpolated g_ij isn't positive definite
+	// CCTK_VWarn() level for interpolated g_ij isn't positive definite
 	// (usually this means we're too close to a singularity)
 	// ... warning level if error occurs on first Newton iteration,
 	//     i.e. when evaluating expansion/Jacobian for initial guess
@@ -177,30 +207,24 @@ struct	error_info
 //
 
 // expansion.cc
-// ... returns true for successful computation,
-// ... returns true for dummy computation (ps_ptr == NULL)
-// ... returns false for failure (eg horizon outside grid);
-//     if (print_msg_flag) then also reports CCTK_VWarn() at the
-//     appropriate warning level
-bool expansion(patch_system* ps_ptr,
-	       const struct cactus_grid_info& cgi,
-	       const struct geometry_info& gi,
-	       const struct error_info& error_info, bool initial_flag,
-	       bool Jacobian_flag = false,
-	       bool print_msg_flag = false,
-	       jtutil::norm<fp>* H_norms_ptr = NULL);
+enum expansion_status
+  expansion(patch_system* ps_ptr,
+	    const struct cactus_grid_info& cgi,
+	    const struct geometry_info& gi,
+	    const struct error_info& error_info, bool initial_flag,
+	    bool Jacobian_flag = false,
+	    bool print_msg_flag = false,
+	    jtutil::norm<fp>* H_norms_ptr = NULL);
 
 // expansion_Jacobian.cc
-// ... returns true for successful computation
-// ... returns true for dummy computation (ps_ptr == NULL && Jac_ptr == NULL)
-// ... returns false for failure
-bool expansion_Jacobian(patch_system* ps_ptr,
-			Jacobian* Jac_ptr,
-			const struct cactus_grid_info& cgi,
-			const struct geometry_info& gi,
-			const struct Jacobian_info& Jacobian_info,
-			const struct error_info& error_info, bool initial_flag,
-			bool print_msg_flag = false);
+enum expansion_status
+  expansion_Jacobian(patch_system* ps_ptr,
+		     Jacobian* Jac_ptr,
+		     const struct cactus_grid_info& cgi,
+		     const struct geometry_info& gi,
+		     const struct Jacobian_info& Jacobian_info,
+		     const struct error_info& error_info, bool initial_flag,
+		     bool print_msg_flag = false);
 
 // Schwarzschild_EF.cc
 void Schwarzschild_EF_geometry(patch_system& ps,