* 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

4 files changed, 262 insertions, 163 deletions

diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index 0ca256e..490ff99 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -56,7 +56,7 @@ using jtutil::pow4;
 
 namespace {
 void setup_xyz_posns(patch_system& ps, bool print_msg_flag);
-bool interpolate_geometry(patch_system& ps,
+bool interpolate_geometry(patch_system* ps_ptr,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag);
@@ -81,9 +81,12 @@ bool compute_Theta(patch_system& ps,
 //******************************************************************************
 
 //
-// This function computes the LHS function Theta(h), and optionally also
-// its Jacobian coefficients (from which the Jacobian matrix may be
-// computed later).
+// If ps_ptr != NULL, this function computes the LHS function Theta(h),
+// and optionally also its Jacobian coefficients (from which the Jacobian
+// matrix may be computed later).
+//
+// If ps_ptr == NULL, this function does a dummy computation, described
+// below.
 //
 // Inputs (angular gridfns, on ghosted grid):
 // ... defined on ghosted grid
@@ -114,6 +117,11 @@ bool compute_Theta(patch_system& ps,
 //	Theta						# $\Theta = \Theta(h)$
 //
 // Arguments:
+// ps_ptr --> The patch system, or == NULL to do (only) a dummy computation
+//	      in which only the parameter-table setup and a dummy geometry
+//	      interpolator call are done, the latter with the number of
+//	      interpolation points is set to 0 and all the output array
+//	      pointers set to NULL.
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
 // print_msg_flag = true to print status messages,
@@ -124,64 +132,88 @@ bool compute_Theta(patch_system& ps,
 //			   to compute various (gridwise) norms of Theta.
 //
 // Results:
-// This function returns true for a successful computation, or 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
+// * 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)
 //
-bool expansion(patch_system& ps,
+bool expansion(patch_system* ps_ptr,
 	       const struct cactus_grid_info& cgi,
 	       const struct geometry_info& gi,
 	       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)
-   then CCTK_VInfo(CCTK_THORNSTRING, "   expansion");
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   %sexpansion",
+		   active_flag ? "" : "dummy ");
 
-// fill in values of all ghosted gridfns in ghost zones
-ps.synchronize();
+if (active_flag)
+   then {
+	//
+	// normal computation
+	//
+
+	// fill in values of all ghosted gridfns in ghost zones
+	ps_ptr->synchronize();
 
-if (gi.check_that_h_is_finite && !h_is_finite(ps, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+	if (gi.check_that_h_is_finite
+	    && !h_is_finite(*ps_ptr, print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
 
-// set up xyz positions of grid points
-setup_xyz_posns(ps, print_msg_flag);
+	// set up xyz positions of grid points
+	setup_xyz_posns(*ps_ptr, print_msg_flag);
+	}
 
 // compute the "geometry" g_ij, K_ij, and partial_k g_ij
 switch	(gi.geometry_method)
 	{
 case geometry__global_interp_from_Cactus_grid:
 case geometry__local_interp_from_Cactus_grid:
-	if (!interpolate_geometry(ps, cgi, gi, print_msg_flag))
+	// 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,
+				  print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
-	if (cgi.Cactus_conformal_metric)
-	   then convert_conformal_to_physical(ps, print_msg_flag);
+	if (active_flag && cgi.Cactus_conformal_metric)
+	   then convert_conformal_to_physical(*ps_ptr, print_msg_flag);
 	break;
 
 case geometry__Schwarzschild_EF:
-	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
+	if (active_flag)
+	   then Schwarzschild_EF_geometry(*ps_ptr, gi, print_msg_flag);
 	break;
+
 default:
 	error_exit(PANIC_EXIT,
 "***** expansion(): unknown gi.geometry_method=(int)%d!\n"
-"                          (this should never happen!)\n"
-,
+"                   (this should never happen!)\n"
+		   ,
 		   int(gi.geometry_method));			/*NOTREACHED*/
 	}
 
-if (gi.check_that_geometry_is_finite
-    && !geometry_is_finite(ps, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+if (active_flag)
+   then {
+	if (gi.check_that_geometry_is_finite
+	    && !geometry_is_finite(*ps_ptr, print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
 
-// compute remaining gridfns --> $\Theta$
-// and optionally also the Jacobian coefficients
-// by algebraic ops and angular finite differencing
-if (!compute_Theta(ps, Jacobian_flag, Theta_norms_ptr, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+	// compute remaining gridfns --> $\Theta$
+	// and optionally also the Jacobian coefficients
+	// by algebraic ops and angular finite differencing
+	if (!compute_Theta(*ps_ptr,
+			   Jacobian_flag, Theta_norms_ptr,
+			   print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
+	}
 
 return true;						// *** NORMAL RETURN ***
 }
@@ -234,7 +266,7 @@ if (print_msg_flag)
 //******************************************************************************
 
 //
-// This function interpolates the Cactus gridfns
+// If ps_ptr != NULL, this function interpolates the Cactus gridfns
 //	gxx...gzz
 //	kxx...kzz
 //	psi			# optional
@@ -245,13 +277,18 @@ if (print_msg_flag)
 //	psi			# optional
 //	partial_d_psi_k		# optional
 // at the nominal-grid trial horizon surface positions given by the
-// global_(x,y,z) angular gridfns.  The psi interpolation is done if and
-// only if the cgi.Cactus_conformal_metric flag is set.  Note that this
-// function ignores the physical-vs-conformal semantics of the gridfns;
-// it just interpolates and takes derivatives of the stored gridfn values.
+// global_(x,y,z) angular gridfns in the patch system *ps_ptr.  The psi
+// interpolation is only done if the cgi.Cactus_conformal_metric flag
+// is set.  Note that this function ignores the physical-vs-conformal
+// semantics of the gridfns; it just interpolates and takes derivatives
+// of the stored gridfn values.
+//
+// If ps_ptr == NULL, this function does (only) the parameter-table
+// setup and a a dummy interpolator call, as described in the comments
+// to  expansion()  above.
 //
-// The interpolation may be either via  CCTK_InterpGridArrays() or via
-// CCTK_InterpLocalUniform() , according to gi.geometry_method.
+// The interpolation may be either via  CCTK_InterpGridArrays()  or via
+//  CCTK_InterpLocalUniform() , according to gi.geometry_method.
 //
 // Inputs (angular gridfns, all on the nominal grid):
 //	global_[xyz]			# xyz positions of grid points
@@ -272,18 +309,23 @@ if (print_msg_flag)
 // This function may also modify the interpolator parameter table.
 //
 // Results:
-// This function returns true for a successful computation, or false
-// if the computation failed because the geometry interpolation would
-// need data outside the Cactus grid, or data from an excised region.
-//
-// FIXME: excision isn't implemented yet :(
+// 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.
 //
 namespace {
-bool interpolate_geometry(patch_system& ps,
+bool interpolate_geometry(patch_system* ps_ptr,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag)
 {
+const bool active_flag = (ps_ptr != NULL);
+
 //
 // Implementation Notes:
 //
@@ -293,8 +335,8 @@ bool interpolate_geometry(patch_system& ps,
 // to either include or exclude the psi entries as appropriate.
 //
 // We remember whether or not psi was interpolated on the previous call,
-// and only modify the interpolator parameter table if this changes (and
-// on our first call).
+// and only modify the interpolator parameter table if this changes (or
+// if this is our first call).
 //
 
 const char* const global_local_str
@@ -311,17 +353,20 @@ if (print_msg_flag)
 
 int status;
 
+#define CAST_PTR_OR_NULL(type_,ptr_)	\
+	(ps_ptr == NULL) ? NULL : static_cast<type_>(ptr_)
+
 
 //
 // ***** interpolation points *****
 //
-const int N_interp_points = ps.N_grid_points();
+const int N_interp_points = (ps_ptr == NULL) ? 0 : ps_ptr->N_grid_points();
 const int interp_coords_type_code = CCTK_VARIABLE_REAL;
 const void* const interp_coords[N_GRID_DIMS]
   = {
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_x)),
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_y)),
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_z)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_x)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_y)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_z)),
     };
 
 
@@ -424,45 +469,45 @@ const CCTK_INT operation_codes[]
 
 void* const output_arrays[]
   = {
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_211)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_311)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_12)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_112)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_212)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_312)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_13)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_113)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_213)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_313)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_22)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_122)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_222)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_322)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_23)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_123)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_223)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_323)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_33)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_133)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_233)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_333)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_11)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_12)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_13)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_33)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__psi)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_1)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_2)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_3)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_11)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_111)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_211)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_311)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_12)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_112)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_212)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_312)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_13)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_113)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_213)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_313)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_22)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_122)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_222)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_322)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_23)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_123)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_223)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_323)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_33)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_133)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_233)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_333)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_11)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_12)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_13)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_22)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_23)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_33)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__psi)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_1)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_2)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_3)),
     };
 
 const int N_output_arrays_for_psi = 4;
-const int N_output_arrays_dim = sizeof(output_arrays)
-				/ sizeof(output_arrays[0]);
+const int N_output_arrays_dim
+	= sizeof(output_arrays) / sizeof(output_arrays[0]);
 const int N_output_arrays_use
 	= psi_flag ? N_output_arrays_dim
 		   : N_output_arrays_dim - N_output_arrays_for_psi;
@@ -526,15 +571,12 @@ if (par_table_setup && (psi_flag == par_table_psi_flag))
 //
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling %s interpolator (%d points)",
-		   global_local_str, N_interp_points);
+		   "         calling %s interpolator (%s%d points)",
+		   global_local_str,
+		   (active_flag ? "" : "dummy: "), N_interp_points);
 switch	(gi.geometry_method)
 	{
 case geometry__global_interp_from_Cactus_grid:
-	if (print_msg_flag)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling CCTK_InterpGridArrays() (%d points)",
-			   N_interp_points);
 	status = CCTK_InterpGridArrays(cgi.GH, N_GRID_DIMS,
 				       gi.operator_handle,
 				       gi.param_table_handle,
@@ -549,10 +591,6 @@ case geometry__global_interp_from_Cactus_grid:
 					  output_arrays);
 	break;
 case geometry__local_interp_from_Cactus_grid:
-	if (print_msg_flag)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling CCTK_InterpLocalUniform() (%d points)",
-			   N_interp_points);
 	status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 					 gi.operator_handle,
 					 gi.param_table_handle,
@@ -581,39 +619,39 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE)
    then {
 	// look in interpolator output table entries
 	// to see *which* point is out-of-range
-	CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end;
+	CCTK_INT error_pt, error_axis, error_direction;
 	if (    (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_pt,
-				  "out_of_range_pt") < 0)
+				  &error_pt,
+				  "error_pt") < 0)
 	     || (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_axis,
-				  "out_of_range_axis") < 0)
+				  &error_axis,
+				  "error_axis") < 0)
 	     || (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_end,
-				  "out_of_range_end") < 0)    )
+				  &error_direction,
+				  "error_direction") < 0)    )
 	   then error_exit(ERROR_EXIT,
 "***** interpolate_geometry():\n"
 "        point out of range when interpolating geometry info from 3-D grid!\n"
 "        (one or more points on the trial horizon surface\n"
 "         are outside the grid or too close to the boundary),\n"
-"        but we can't get info about the out-of-range point!\n"
+"        but we can't get info about the bad point(s)!\n"
 "        ==> maybe an interpolator problem?\n");		/*NOTREACHED*/
 
-	assert(out_of_range_pt >= 0);
-	assert(out_of_range_pt < ps.N_grid_points());
+	assert(error_pt >= 0);
+	assert(error_pt < ps_ptr->N_grid_points());
 	const double global_x
-	   = ps.gridfn_data(gfns::gfn__global_x)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_x)[error_pt];
 	const double global_y
-	   = ps.gridfn_data(gfns::gfn__global_y)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_y)[error_pt];
 	const double global_z
-	   = ps.gridfn_data(gfns::gfn__global_z)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_z)[error_pt];
 
-	assert(out_of_range_axis >= 0);
-	assert(out_of_range_axis < N_GRID_DIMS);
-	const char axis = "xyz"[out_of_range_axis];
+	assert(error_axis >= 0);
+	assert(error_axis < N_GRID_DIMS);
+	const char axis = "xyz"[error_axis];
 
-	assert((out_of_range_end == -1) || (out_of_range_end == +1));
-	const char end = (out_of_range_end == -1) ? '-' : '+';
+	assert((error_direction == -1) || (error_direction == +1));
+	const char end = (error_direction == -1) ? '-' : '+';
 
 	if (print_msg_flag)
 	   then {
@@ -636,7 +674,6 @@ if (status < 0)
 ,
 		   status);					/*NOTREACHED*/
 
-
 return true;						// *** NORMAL RETURN ***
 }
 	  }
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index 082387b..805ed97 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -67,8 +67,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(patch_system& ps,
-			      Jacobian& Jac,
+bool expansion_Jacobian_dr_FD(patch_system* ps_ptr,
+			      Jacobian* Jac_ptr,
 			      const struct cactus_grid_info& cgi,
 			      const struct geometry_info& gi,
 			      const struct Jacobian_info& Jacobian_info,
@@ -78,32 +78,57 @@ bool expansion_Jacobian_dr_FD(patch_system& ps,
 //******************************************************************************
 
 //
-// This function is a top-level driver to compute the Jacobian matrix
-// J[Theta(h)] of the expansion Theta(h).  It just decodes the
-// Jacobian_compute_method parameter and calls the appropriate subfunction.
+// If ps_ptr != NULL and Jac_ptr != NULL, this function computes the
+// Jacobian matrix J[Theta(h)] of the expansion Theta(h).  We assume
+// that Theta(h) has already been computed.
 //
-// We assume that Theta(h) has already been computed.
+// If ps_ptr == NULL and Jac_ptr == NULL, this function does a dummy
+// computation, in which only any expansion() (and hence geometry
+// interpolator) calls are done, these with the number of interpolation
+// points set to 0 and all the output array pointers set to NULL.
+//
+// It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL.
+//
+// Only some values of  Jacobian_info.Jacobian_compute_method  support
+// the dummy computation.
+//
+// Arguments:
+// ps_ptr --> The patch system, or == NULL to do (only) a dummy computation.
+// Jac_ptr --> The Jacobian, or == NULL to do (only) a dummy computation.
 //
 // 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
+// * 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,
-			Jacobian& Jac,
+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,
 			bool print_msg_flag /* = false */)
 {
+const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL);
+
 switch	(Jacobian_info.Jacobian_compute_method)
 	{
 case Jacobian__numerical_perturbation:
-	if (!  expansion_Jacobian_NP(ps, Jac,
-				   cgi, gi, Jacobian_info,
-				   print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
-	break;
+	if (active_flag)
+	   then {
+		if (! expansion_Jacobian_NP(*ps_ptr, *Jac_ptr,
+					    cgi, gi, Jacobian_info,
+					    print_msg_flag))
+		   then return false;			// *** ERROR RETURN ***
+		break;
+		}
+	   else CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   expansion_Jacobian():\n"
+"        dummy computation isn't supported for\n"
+"        Jacobian_compute_method = \"numerical perturbation\"!");/*NOTREACHED*/
+
 case Jacobian__symbolic_diff:
 	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -112,15 +137,20 @@ case Jacobian__symbolic_diff:
 "        isn't implemented (yet); reverting to\n"
 "        \"symbolic differentiation with finite diff d/dr\"");
 	// fall through
+
 case Jacobian__symbolic_diff_with_FD_dr:
-	expansion_Jacobian_partial_SD(ps, Jac,
-				      cgi, gi, Jacobian_info,
-				      print_msg_flag);
-	if (! expansion_Jacobian_dr_FD(ps, Jac,
+	if (active_flag)
+	   then expansion_Jacobian_partial_SD(*ps_ptr, *Jac_ptr,
+					      cgi, gi, Jacobian_info,
+					      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,
 				       print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
 	break;
+
 default:
 	error_exit(PANIC_EXIT,
 "***** expansion_Jacobian():\n"
@@ -129,6 +159,7 @@ default:
 ,
 		   int(Jacobian_info.Jacobian_compute_method));	/*NOTREACHED*/
 	}
+
 return true;						// *** NORMAL RETURN ***
 }
 
@@ -203,8 +234,9 @@ 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, ii))
+	if (! expansion(&ps, cgi, ii))
 	   then return false;				// *** ERROR RETURN ***
+
 		for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
 		{
 		patch& xp = ps.ith_patch(xpn);
@@ -226,6 +258,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 		}
 		}
 		}
+
 	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y;
 	}
 	}
@@ -459,8 +492,16 @@ patch& yp = ye.my_patch();
 //******************************************************************************
 
 //
-// This function sums the d/dr terms into the Jacobian matrix of the
-// expansion Theta(h), computing those terms by finite differencing.
+// If ps_ptr != NULL and Jac_ptr != NULL, this function sums the d/dr
+// terms into the Jacobian matrix of the expansion Theta(h), computing
+// those terms by finite differencing.
+//
+// If ps_ptr == NULL and Jac_ptr == NULL, this function does a dummy
+// computation, in which only any expansion() (and hence geometry
+// interpolator) calls are done, these with the number of interpolation
+// points set to 0 and all the output array pointers set to NULL.
+//
+// It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL.
 //
 // The basic algorithm is that
 //	Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ]
@@ -473,53 +514,65 @@ patch& yp = ye.my_patch();
 //	Jac += d/dr terms
 //
 // Results:
-// This function returns true for a successful computation, or false
-// if the computation failed because the geometry interpolation would
-// need data outside the Cactus grid, or data from an excised region.
-// FIXME: excision isn't implemented yet :(
+// 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 :(
 //
 namespace {
-bool expansion_Jacobian_dr_FD(patch_system& ps,
-			      Jacobian& Jac,
+bool expansion_Jacobian_dr_FD(patch_system* ps_ptr,
+			      Jacobian* Jac_ptr,
 			      const struct cactus_grid_info& cgi,
 			      const struct geometry_info& gi,
 			      const struct Jacobian_info& Jacobian_info,
 			      bool print_msg_flag)
 {
+const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL);
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   horizon Jacobian: d/dr terms (finite diff)");
+		   "   horizon Jacobian: %sd/dr terms (finite diff)",
+		   active_flag ? "" : "dummy ");
 
 const fp epsilon = Jacobian_info.perturbation_amplitude;
 
 // compute Theta(h+epsilon)
-ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
-ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
-if (! expansion(ps, cgi, gi))
+if (active_flag)
+   then {
+	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))
    then return false;					// *** ERROR RETURN ***
 
-	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
-	{
-	patch& p = ps.ith_patch(pn);
+if (active_flag)
+   then {
+	    for (int pn = 0 ; pn < ps_ptr->N_patches() ; ++pn)
+	    {
+	    patch& p = ps_ptr->ith_patch(pn);
 		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
 		{
 		for (int isigma = p.min_isigma() ;
 		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
-		const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma);
+		const int II = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma);
 		const fp old_Theta = p.gridfn(gfns::gfn__save_Theta,
 					      irho,isigma);
 		const fp new_Theta = p.gridfn(gfns::gfn__Theta,
 					      irho,isigma);
-		Jac.sum_into_element(II,II, (new_Theta - old_Theta) / epsilon);
+		const fp d_dr_term = (new_Theta - old_Theta) / epsilon;
+		Jac_ptr->sum_into_element(II,II, d_dr_term);
 		}
 		}
-	}
+	    }
 
-// restore h and Theta
-ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
-ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
+	// restore h and Theta
+	ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+	ps_ptr->gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
+	}
 
 return true;						// *** NORMAL RETURN ***
 }
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index c20a720..3014cd7 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -144,6 +144,7 @@ struct	Jacobian_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
@@ -156,7 +157,7 @@ namespace expansion_warning_levels
 					// ==> skipping nonfinite check
 	       };
 	  }
-bool expansion(patch_system& ps,
+bool expansion(patch_system* ps_ptr,
 	       const struct cactus_grid_info& cgi,
 	       const struct geometry_info& gi,
 	       bool Jacobian_flag = false,
@@ -164,9 +165,11 @@ bool expansion(patch_system& ps,
 	       jtutil::norm<fp>* H_norms_ptr = NULL);
 
 // expansion_Jacobian.cc
-// ... returns true for successful computation, false for failure
-bool expansion_Jacobian(patch_system& ps,
-			Jacobian& Jac,
+// ... 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,
@@ -180,6 +183,8 @@ void Schwarzschild_EF_geometry(patch_system& ps,
 // misc.cc
 enum geometry_method
   decode_geometry_method(const char geometry_method_string[]);
+#ifdef NOT_USED
 bool geometry_method_is_interp(enum geometry_method geometry_method);
+#endif
 enum Jacobian_compute_method
   decode_Jacobian_compute_method(const char Jacobian_compute_method_string[]);
diff --git a/src/gr/misc.cc b/src/gr/misc.cc
index 7cdc018..8a0ced4 100644
--- a/src/gr/misc.cc
+++ b/src/gr/misc.cc
@@ -2,7 +2,9 @@
 // $Header$
 //
 // decode_geometry_method - decode the geometry_method parameter
+#ifdef NOT_USED
 // geomery_method_is_interp - does  enum geometry_method  specify interpolation?
+#endif
 // decode_geometry_method - decode the geometry_method parameter
 //
 
@@ -61,6 +63,7 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 
 //******************************************************************************
 
+#ifdef NOT_USED
 //
 // This function determines whether or not an  enum geometry_method
 // value specifies an interpolation.
@@ -80,6 +83,7 @@ default:
 		   geometry_method);				/*NOTREACHED*/
 	}
 }
+#endif
 
 //******************************************************************************