* add support for CCTK_InterpGridArrays() returning per-point

  interpolator status (use it if it's available, otherwise fall back
  to CCTK_InterpGridArrays() function result)
* drop support for interpolators returning status info on which side
  of the grid an outside-the-grid point is outside -- this turned out
  to be awkward to implement in combinatino with per-point status,
  I only used it to produce slightly more informative error messages
  in a few cases (ditto in ../patch/)
* convert some commented-out debugging code to #ifdef on various DEBUG symbols


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

1 files changed, 139 insertions, 64 deletions

diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index f26cac1..d6c51b5 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -14,6 +14,26 @@
 /// compute_Theta - compute Theta(h) given earlier setup
 ///
 
+//
+// debugging flags
+//
+#undef GEOMETRY_INTERP_DEBUG	// define this for verbose debugging
+				// of geometry interpolator calls/results
+#undef GEOMETRY_INTERP_DEBUG2	// define this for even more verbose debugging
+				// of geometry interpolator calls/results
+#undef GEOMETRY_INTERP_SYNC_SLEEP	// use sleep() to try to ensure we're
+					// synchronized across all processors
+					// before calling geometry interpolator
+#undef GEOMETRY_INTERP_SYNC_BARRIER	// use CCTK_Barrier() to ensure we're
+					// synchronized across all processors
+					// before calling geometry interpolator
+#undef COMPUTE_THETA_DEBUG	// define this for verbose debugging
+				// in compute_Theta()
+
+#ifdef GEOMETRY_INTERP_SYNC_SLEEP
+#include <unistd.h>		// sleep()
+#endif
+
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
@@ -190,12 +210,22 @@ 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
 	  {
+	#ifdef GEOMETRY_INTERP_DEBUG
+	printf("AHFinderDirect:: proc %d: about to call interpolate_geometry()\n",
+	       int(CCTK_MyProc(cgi.GH)));
+	fflush(stdout);
+	#endif
 	const
 	  enum expansion_status status
 	    = interpolate_geometry(ps_ptr,
 				   cgi, gi,
 				   error_info, initial_flag,
 				   print_msg_flag);
+	#ifdef GEOMETRY_INTERP_DEBUG
+	printf("AHFinderDirect:: proc %d: interpolate_geometry() status=(int)%d\n",
+	       int(CCTK_MyProc(cgi.GH)), int(status));
+	fflush(stdout);
+	#endif
 	if (status != expansion_success)
 	   then return status;				// *** ERROR RETURN ***
 	  }
@@ -614,9 +644,38 @@ if (print_msg_flag)
 		   "         calling %s interpolator (%s%d points)",
 		   global_local_str,
 		   (active_flag ? "" : "dummy: "), N_interp_points);
+
+#ifdef GEOMETRY_INTERP_DEBUG
+printf("AHFinderDirect:: proc %d: initializing interpolator outputs to 999.999\n",
+       int(CCTK_MyProc(cgi.GH)));
+  {
+	for (int pt = 0 ; pt < N_interp_points ; ++pt)
+	{
+		for (int out = 0 ; out < N_output_arrays_use ; ++out)
+		{
+		CCTK_REAL* const out_ptr
+			= static_cast<CCTK_REAL*>(output_arrays[out]);
+		out_ptr[pt] = 999.999;
+		}
+	}
+  }
+#endif
+
 switch	(gi.geometry_method)
 	{
 case geometry__global_interp_from_Cactus_grid:
+	  {
+	#ifdef GEOMETRY_INTERP_DEBUG
+	printf("AHFinderDirect:: proc %d: calling global interpolator (N_interp_points=%d)\n",
+	       int(CCTK_MyProc(cgi.GH)), N_interp_points);
+	fflush(stdout);
+	#endif
+	#ifdef GEOMETRY_INTERP_SYNC_SLEEP
+	sleep(1);
+	#endif
+	#ifdef GEOMETRY_INTERP_SYNC_BARRIER
+	CCTK_Barrier(cgi.GH);
+	#endif
 	status = CCTK_InterpGridArrays(cgi.GH, N_GRID_DIMS,
 				       gi.operator_handle,
 				       gi.param_table_handle,
@@ -629,7 +688,64 @@ case geometry__global_interp_from_Cactus_grid:
 				       N_output_arrays_use,
 					  output_array_type_codes,
 					  output_arrays);
+	#ifdef GEOMETRY_INTERP_DEBUG
+	printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() returned status=%d\n",
+	       int(CCTK_MyProc(cgi.GH)), status);
+	fflush(stdout);
+	#endif
+	#ifdef GEOMETRY_INTERP_DEBUG2
+	  {
+		for (int pt = 0 ; pt < N_interp_points ; pt = 2*pt + (pt == 0))
+		{
+		printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() results for pt=%d:\n",
+		       int(CCTK_MyProc(cgi.GH)), pt);
+			for (int out = 0 ; out < N_output_arrays_use ; ++out)
+			{
+			const CCTK_REAL* const out_ptr
+				= static_cast<const CCTK_REAL*>(output_arrays[out]);
+			printf("   out=%d   result=%g\n", out, double(out_ptr[pt]));
+			}
+		}
+	  }
+	#endif	/* GEOMETRY_INTERP_DEBUG2 */
+	//
+	// CCTK_InterpGridArrays() returns a status that reflects
+	// the interpolation on *all* processors, but what we really
+	// want here is just the status for the interpolation of the
+	// points that *we* (this processor) asked for ==> get that
+	// status (assuming the global and local interpolators support
+	// providing it)
+	//
+	CCTK_INT local_interpolator_status;
+	const int table_status = Util_TableGetInt(gi.param_table_handle,
+						  &local_interpolator_status,
+						  "local_interpolator_status");
+	if	(table_status == 1)
+	   then {
+		// we got the local interpolator status successfully
+		// ==> use it for all further checks
+		#ifdef GEOMETRY_INTERP_DEBUG
+		printf("AHFinderDirect:: proc %d: replacing status=%d with local_interpolator_status=%d\n",
+		       int(CCTK_MyProc(cgi.GH)), status, int(local_interpolator_status));
+		fflush(stdout);
+		#endif
+		status = local_interpolator_status;
+		}
+	else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+	   then {
+		// evidently the interpolators don't provide the
+		// local interpolator status
+		// ==> stick with the status have
+		// ==> no-op here
+		}
+	else	error_exit(ERROR_EXIT,
+"***** interpolate_geometry():\n"
+"        error return %d trying to get local interpolator status\n"
+"        from parameter table!  (CCTK_InterpGridArrays() status=%d)\n"
+			   ,
+			   table_status, status);		/*NOTREACHED*/
 	break;
+	  }
 case geometry__local_interp_from_Cactus_grid:
 	status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 					 gi.operator_handle,
@@ -660,77 +776,17 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE)
    then {
 	if (print_msg_flag)
 	   then {
-		// did the local interpolator report an error on this processor?
+		// see if we can get further info
 		const int warn_level
 		   = initial_flag
 		     ? error_info.warn_level__point_outside__initial
 		     : error_info.warn_level__point_outside__subsequent;
 
-		// look in interpolator output table entries
-		// to see if the local interpolator reported an error on this
-		// processor, and if so, *which* point is out-of-range
-		CCTK_INT local_interpolator_status;
-		const int status_interp
-		   = Util_TableGetInt(gi.param_table_handle,
-				      &local_interpolator_status,
-				      "local_interpolator_status");
-		const bool point_outside_on_this_processor
-		   = (status_interp > 0)
-		     && (local_interpolator_status == CCTK_ERROR_INTERP_POINT_OUTSIDE);
-
-		CCTK_INT error_pt, error_axis, error_direction;
-		const int status_pt
-		   = Util_TableGetInt(gi.param_table_handle,
-				      &error_pt, "error_pt");
-		const int status_axis
-		   = Util_TableGetInt(gi.param_table_handle,
-				      &error_axis, "error_axis");
-		const int status_direction
-		   = Util_TableGetInt(gi.param_table_handle,
-				      &error_direction, "error_direction");
-		if (point_outside_on_this_processor && (status_pt > 0)
-						    && (status_axis > 0)
-						    && (status_direction > 0))
-		   then {
-			assert(error_pt >= 0);
-			assert(error_pt < ps_ptr->N_grid_points());
-			const double global_x
-			   = ps_ptr->gridfn_data(gfns::gfn__global_x)[error_pt];
-			const double global_y
-			   = ps_ptr->gridfn_data(gfns::gfn__global_y)[error_pt];
-			const double global_z
-			   = ps_ptr->gridfn_data(gfns::gfn__global_z)[error_pt];
-
-			assert(error_axis >= 0);
-			assert(error_axis < N_GRID_DIMS);
-			const char axis = "xyz"[error_axis];
-
-			assert(    (error_direction == -1)
-				|| (error_direction == +1)    );
-			const char direction
-			   = (error_direction == -1) ? '-' : '+';
-
-			CCTK_VWarn(warn_level, __LINE__, __FILE__,
-				   CCTK_THORNSTRING,
-"\n"
-"interpolate_geometry():\n"
-"        the trial-horizon-surface point [pt=%d]\n"
-"           global_(x,y,z)=(%g,%g,%g)\n"
-"        is outside the grid (or too close to the boundary)"
-						" in the %c%c direction!\n"
-"        (in a single-processor run, this may also mean that\n"
-"         driver::ghost_size is too small for this geometry interpolator)\n"
-				   ,
-				   error_pt,
-				   global_x, global_y, global_z,
-				   direction, axis);
-			}
-		   else CCTK_VWarn(warn_level, __LINE__, __FILE__,
-				   CCTK_THORNSTRING,
-"\n"
+
+		CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING,
 "interpolate_geometry():\n"
 "        one or more points on the trial horizon surface point"
-"        is/are outside the grid (or too close to the boundary)"
+"        is/are outside the grid (or too close to the grid boundary)"
 "        (in a single-processor run, this may also mean that\n"
 "         driver::ghost_size is too small for this geometry interpolator)\n");
 		}
@@ -1179,6 +1235,25 @@ if (print_msg_flag)
 		const fp sigma = p.sigma_of_isigma(isigma);
 		fp xx, yy, zz;
 		p.xyz_of_r_rho_sigma(r,rho,sigma, xx, yy, zz);
+		#ifdef COMPUTE_THETA_DEBUG
+		if (    ((irho == 0) && (isigma == 0))
+		     || ((irho == 5) && (isigma == 5))    )
+		   then {
+			printf("AHFinderDirect:: computing at %s patch (%d,%d) r=%g ==> local xyz=(%g,%g,%g)\n",
+			       p.name(), irho,isigma, double(r),
+			       double(xx), double(yy), double(zz));
+			printf("AHFinderDirect:: got g_dd_11=%g g_dd_12=%g g_dd_13=%g\n",
+			       p.gridfn(gfns::gfn__g_dd_11, irho,isigma),
+			       p.gridfn(gfns::gfn__g_dd_12, irho,isigma),
+			       p.gridfn(gfns::gfn__g_dd_13, irho,isigma));
+			printf("                                g_dd_22=%g g_dd_23=%g\n",
+			       p.gridfn(gfns::gfn__g_dd_22, irho,isigma),
+			       p.gridfn(gfns::gfn__g_dd_23, irho,isigma));
+			printf("                                           g_dd_33=%g\n",
+			       p.gridfn(gfns::gfn__g_dd_33, irho,isigma));
+			fflush(stdout);
+			}
+		#endif
 
 		// 1st derivative coefficients X_ud
 		const fp X_ud_11 = p.partial_rho_wrt_x(xx, yy, zz);