add beginnings of support for using new CCTK_InterpGridArrays()

global interpolator & multiprocessor operation

... at present the CCTK_InterpGridArrays() support works fine,
    but only on 1 processor, I seem to get an infinite loop on
    multiple processors


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

10 files changed, 209 insertions, 135 deletions

diff --git a/par/Kerr.par b/par/Kerr.par
index d57121d..003c5ad 100644
--- a/par/Kerr.par
+++ b/par/Kerr.par
@@ -3,24 +3,23 @@
 # and the initial guess are both deliberately de-centered with respect
 # to the black hole, to make this a non-trivial test for the apparent
 # horizon finder.
-ActiveThorns = "CartGrid3D LocalInterp PUGH ADMBase ADMCoupling StaticConformal CoordGauge Exact AHFinderDirect"
 
 # flesh
 cactus::cctk_itlast = 0
 
-# PUGH
+ActiveThorns = "PUGH"
 driver::ghost_size = 2
 driver::global_nx = 31
 driver::global_ny = 31
 driver::global_nz = 17
 
-# CartGrid3D
+ActiveThorns = "CartGrid3D"
 grid::domain = "bitant"
 grid::avoid_origin = "false"
 grid::type = "byspacing"
 grid::dxyz = 0.2
 
-# ADMBase
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
 ADMBase::initial_lapse = "exact"
 ADMBase::initial_shift = "exact"
 ADMBase::initial_data = "exact"
@@ -28,13 +27,14 @@ ADMBase::lapse_evolution_method = "static"
 ADMBase::shift_evolution_method = "static"
 ADMBase::metric_type = "physical"
 
-# Exact
 Exact::exact_model = "Kerr/Kerr-Schild"
 Exact::Kerr_KerrSchild__mass = 1.0
 Exact::Kerr_KerrSchild__spin = 0.6
 
 ########################################
 
+ActiveThorns = "LocalInterp PUGHInterp AHFinderDirect"
+
 AHFinderDirect::find_AHs_at_postinitial = "true"
 AHFinderDirect::find_AHs_at_poststep    = "false"
 AHFinderDirect::print_timing_stats = "true"
diff --git a/param.ccl b/param.ccl
index 9b8ebe5..8466572 100644
--- a/param.ccl
+++ b/param.ccl
@@ -545,13 +545,20 @@ keyword integral_method "how do we compute integrals over the horizon?"
 
 keyword geometry_method "how do we compute the slice's geometry?"
 {
-# for normal use
-"interpolate from Cactus grid" :: "interpolate gij and Kij from Cactus grid"
+# this would be for normal use, but it doesn't work yet :( :( :(
+"globally interpolate from Cactus grid" :: \
+   "interpolate gij and Kij from global Cactus grid \
+    using CCTK_InterpGridArrays() global interpolator API"
+
+# use this for now
+"locally interpolate from Cactus grid" :: \
+   "interpolate gij and Kij from local Cactus grid on this processor \
+    using CCTK_InterpLocalUniform() localinterpolator API"
 
 # for testing/debugging
 "Schwarzschild/EF"             :: \
   "hard-wire to Schwarzschild spacetime / Eddington-Finkelstein slice"
-} "interpolate from Cactus grid"
+} "locally interpolate from Cactus grid"
 
 ########################################
 
@@ -562,10 +569,6 @@ keyword geometry_method "how do we compute the slice's geometry?"
 # Cactus grid to the position of each trial horizon surface, giving
 # gij, Kij, and partial_x gij as outputs.  This interpolator must have
 # the following properties:
-# - It must be able to accept the Cactus grid as an input domain.  The
-#   current implementation uses the  CCTK_InterpLocalUniform()  API; this
-#   is processor-local and so restricts the whole apparent horizon finder
-#   to single-processor operation.
 # - It must support taking at least 1st derivatives as part of the
 #   interpolation.
 # - It should give at least $C^1$ interpolants for smooth data, otherwise
@@ -578,10 +581,16 @@ keyword geometry_method "how do we compute the slice's geometry?"
 # thorn).
 #
 
+string coordinate_system_name \
+  "name under which the coordinate system is registered in Cactus"
+{
+.+ :: "any string (in practice it should be nonempty)"
+} "cart3d"
+
 string geometry_interpolator_name \
   "name under which the geometry interpolation operator is registered in Cactus"
 {
-.* :: "any string"
+.+ :: "any string"
 } "Hermite polynomial interpolation"
 
 string geometry_interpolator_pars \
@@ -679,7 +688,7 @@ boolean check_that_geometry_is_finite \
 string interpatch_interpolator_name \
   "name under which the interpatch interpolation operator is registered in Cactus"
 {
-.* :: "any string"
+.+ :: "any string (in practice it should be nonempty)"
 } "Lagrange polynomial interpolation"
 
 string interpatch_interpolator_pars \
diff --git a/run/test-ahfinderdirect/misc/Kerr.par b/run/test-ahfinderdirect/misc/Kerr.par
index 9ba6be0..003c5ad 100644
--- a/run/test-ahfinderdirect/misc/Kerr.par
+++ b/run/test-ahfinderdirect/misc/Kerr.par
@@ -1,22 +1,25 @@
-# parameter file for patch system test
-ActiveThorns = "CartGrid3D LocalInterp PUGH ADMBase ADMCoupling StaticConformal CoordGauge Exact AHFinderDirect"
+# This parameter file sets up Kerr/Kerr-Schild initial data, then
+# finds the apparent horizon in it.  The local coordinate system origin
+# and the initial guess are both deliberately de-centered with respect
+# to the black hole, to make this a non-trivial test for the apparent
+# horizon finder.
 
 # flesh
 cactus::cctk_itlast = 0
 
-# PUGH
-driver::ghost_size = 0
-driver::global_nx = 61
-driver::global_ny = 61
-driver::global_nz = 61
+ActiveThorns = "PUGH"
+driver::ghost_size = 2
+driver::global_nx = 31
+driver::global_ny = 31
+driver::global_nz = 17
 
-# CartGrid3D
-grid::domain = "full"
-grid::type = "byrange"
-grid::xyzmin = -3.0
-grid::xyzmax =  3.0
+ActiveThorns = "CartGrid3D"
+grid::domain = "bitant"
+grid::avoid_origin = "false"
+grid::type = "byspacing"
+grid::dxyz = 0.2
 
-# ADMBase
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
 ADMBase::initial_lapse = "exact"
 ADMBase::initial_shift = "exact"
 ADMBase::initial_data = "exact"
@@ -24,17 +27,30 @@ ADMBase::lapse_evolution_method = "static"
 ADMBase::shift_evolution_method = "static"
 ADMBase::metric_type = "physical"
 
-# Exact
 Exact::exact_model = "Kerr/Kerr-Schild"
 Exact::Kerr_KerrSchild__mass = 1.0
 Exact::Kerr_KerrSchild__spin = 0.6
 
-AHFinderDirect::find_AHs_at_poststep = "false"
+########################################
 
-AHFinderDirect::patch_system_type = "full sphere"
-AHFinderDirect::N_zones_per_right_angle = 18
+ActiveThorns = "LocalInterp PUGHInterp AHFinderDirect"
+
+AHFinderDirect::find_AHs_at_postinitial = "true"
+AHFinderDirect::find_AHs_at_poststep    = "false"
+AHFinderDirect::print_timing_stats = "true"
+
+AHFinderDirect::final_Theta_update_if_Delta_h_converged = "true"
+AHFinderDirect::how_often_to_output_Theta = 1
+AHFinderDirect::h_base_file_name     = "Kerr.h"
+AHFinderDirect::Theta_base_file_name = "Kerr.Theta"
 
 AHFinderDirect::N_horizons = 1
-AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild"
-AHFinderDirect::initial_guess__Kerr_KerrSchild__mass[1] = 1.0
-AHFinderDirect::initial_guess__Kerr_KerrSchild__spin[1] = 0.6
+AHFinderDirect::origin_x[1] = 0.5
+AHFinderDirect::origin_y[1] = 0.7
+AHFinderDirect::origin_z[1] = 0.0
+
+AHFinderDirect::initial_guess_method = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2
+AHFinderDirect::initial_guess__coord_sphere__y_center[1] =  0.3
+AHFinderDirect::initial_guess__coord_sphere__z_center[1] =  0.0
+AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 87afdf5..9cd0d1f 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -199,6 +199,7 @@ struct	state
 	enum method method;
 	struct verbose_info verbose_info;
 	int timer_handle;
+	int my_proc;			// our processor number
 
 	struct IO_info IO_info;
 	struct Jacobian_info Jac_info;
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index c9ef0b5..e840480 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -146,10 +146,12 @@ const bool output_Theta
      && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0);
 
 //
-// 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 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
 //
-setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
+if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
+   then setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 
 	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
 	{
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index d48fbe7..5151740 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -118,6 +118,7 @@ verbose_info.print_algorithm_details
    = (state.verbose_info.verbose_level >= verbose_level__algorithm_details);
 
 state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1;
+state.my_proc = CCTK_MyProc(cctkGH);
 
 struct IO_info& IO_info = state.IO_info;
 IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
@@ -164,16 +165,13 @@ if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
 struct cactus_grid_info& cgi = state.cgi;
 cgi.GH = cctkGH;
-cgi.coord_origin[0] = cctk_origin_space[0];
-cgi.coord_origin[1] = cctk_origin_space[1];
-cgi.coord_origin[2] = cctk_origin_space[2];
-cgi.coord_delta[0] = cctk_delta_space[0];
-cgi.coord_delta[1] = cctk_delta_space[1];
-cgi.coord_delta[2] = cctk_delta_space[2];
-cgi.gridfn_dims[0] = cctk_lsh[0];
-cgi.gridfn_dims[1] = cctk_lsh[1];
-cgi.gridfn_dims[2] = cctk_lsh[2];
+cgi.Cactus_conformal_metric = false;	// dummy value, may change later
 
+cgi.coord_system_handle = CCTK_CoordSystemHandle(coordinate_system_name);
+if (cgi.coord_system_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_setup(): can't find Cactus coordinate system \"%\"!",
+		   coordinate_system_name);			/*NOTREACHED*/
 cgi.g_dd_11_varindex = CCTK_VarIndex("ADMBase::gxx");
 cgi.g_dd_12_varindex = CCTK_VarIndex("ADMBase::gxy");
 cgi.g_dd_13_varindex = CCTK_VarIndex("ADMBase::gxz");
@@ -188,16 +186,42 @@ cgi.K_dd_23_varindex = CCTK_VarIndex("ADMBase::kyz");
 cgi.K_dd_33_varindex = CCTK_VarIndex("ADMBase::kzz");
 cgi.psi_varindex     = CCTK_VarIndex("StaticConformal::psi");
 
+cgi.coord_origin[0] = cctk_origin_space[0];
+cgi.coord_origin[1] = cctk_origin_space[1];
+cgi.coord_origin[2] = cctk_origin_space[2];
+cgi.coord_delta[0] = cctk_delta_space[0];
+cgi.coord_delta[1] = cctk_delta_space[1];
+cgi.coord_delta[2] = cctk_delta_space[2];
+cgi.gridfn_dims[0] = cctk_lsh[0];
+cgi.gridfn_dims[1] = cctk_lsh[1];
+cgi.gridfn_dims[2] = cctk_lsh[2];
+
+// we can't set the data pointers yet because they may change from time
+// step to time step as (if) Cactus rotates the multiple time levels
+cgi.g_dd_11_data = NULL;		// dummy value, will be reset later
+cgi.g_dd_12_data = NULL;		// dummy value, will be reset later
+cgi.g_dd_13_data = NULL;		// dummy value, will be reset later
+cgi.g_dd_22_data = NULL;		// dummy value, will be reset later
+cgi.g_dd_23_data = NULL;		// dummy value, will be reset later
+cgi.g_dd_33_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_11_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_12_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_13_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_22_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_23_data = NULL;		// dummy value, will be reset later
+cgi.K_dd_33_data = NULL;		// dummy value, will be reset later
+cgi.psi_data     = NULL;		// dummy value, will be reset later
+
 
 //
 // set up the geometry parameters and the geometry interpolator
 //
+if (verbose_info.print_algorithm_highlights)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
 struct geometry_info& gi = state.gi;
 gi.geometry_method = decode_geometry_method(geometry_method);
 
-// parameters for geometry_method = "interpolate from Cactus grid"
-if (verbose_info.print_algorithm_highlights)
-   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
+// parameters for geometry_method = interpolation methods
 gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
 if (gi.operator_handle < 0)
    then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index 9cafc2b..f2b5961 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -3,7 +3,6 @@
 //
 // <<<prototypes for functions local to this file>>>
 // expansion - top-level driver
-// decode_geometry_method - decode the geometry_method parameter
 //
 /// setup_xyz_posns - setup global xyz posns of grid points
 /// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid
@@ -148,8 +147,7 @@ if (print_msg_flag)
 // fill in values of all ghosted gridfns in ghost zones
 ps.synchronize();
 
-if (gi.check_that_h_is_finite
-    && !h_is_finite(ps, print_msg_flag))
+if (gi.check_that_h_is_finite && !h_is_finite(ps, print_msg_flag))
    then return false;					// *** ERROR RETURN ***
 
 // set up xyz positions of grid points
@@ -158,12 +156,14 @@ setup_xyz_posns(ps, print_msg_flag);
 // compute the "geometry" g_ij, K_ij, and partial_k g_ij
 switch	(gi.geometry_method)
 	{
-case geometry__interpolate_from_Cactus_grid:
+case geometry__global_interp_from_Cactus_grid:
+case geometry__local_interp_from_Cactus_grid:
 	if (!interpolate_geometry(ps, cgi, gi, print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
 	if (cgi.Cactus_conformal_metric)
 	   then convert_conformal_to_physical(ps, print_msg_flag);
 	break;
+
 case geometry__Schwarzschild_EF:
 	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
 	break;
@@ -189,24 +189,6 @@ return true;						// *** NORMAL RETURN ***
 }
 
 //******************************************************************************
-
-//
-// This function decodes the  geometry_method  parameter (string) into
-// an internal enum for future use.
-//
-enum geometry_method
-  decode_geometry_method(const char geometry_method_string[])
-{
-if	(STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid"))
-   then return geometry__interpolate_from_Cactus_grid;
-else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF"))
-   then return geometry__Schwarzschild_EF;
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"decode_geometry_method(): unknown geometry_method_string=\"%s\"!",
-		   geometry_method_string);			/*NOTREACHED*/
-}
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -270,6 +252,9 @@ if (print_msg_flag)
 // function ignores the physical-vs-conformal semantics of the gridfns;
 // it just interpolates and takes derivatives of the stored gridfn values.
 //
+// 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
 //
@@ -314,12 +299,17 @@ bool interpolate_geometry(patch_system& ps,
 // on our first call).
 //
 
+const char* const global_local_str
+	= (gi.geometry_method == geometry__global_interp_from_Cactus_grid)
+	  ? "global"
+	  : "local";
 const bool psi_flag = cgi.Cactus_conformal_metric;
 
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      interpolating %s from Cactus grid",
-		   psi_flag ? "{g_ij, K_ij, psi}" : "{g_ij, K_ij}");
+		   global_local_str,
+		   (psi_flag ? "{g_ij, K_ij, psi}" : "{g_ij, K_ij}"));
 
 int status;
 
@@ -341,6 +331,16 @@ const void* const interp_coords[N_GRID_DIMS]
 // ***** input arrays *****
 //
 
+const CCTK_INT input_array_variable_indices[]
+	= {
+	  cgi.g_dd_11_varindex, cgi.g_dd_12_varindex, cgi.g_dd_13_varindex,
+				cgi.g_dd_22_varindex, cgi.g_dd_23_varindex,
+						      cgi.g_dd_33_varindex,
+	  cgi.K_dd_11_varindex, cgi.K_dd_12_varindex, cgi.K_dd_13_varindex,
+				cgi.K_dd_22_varindex, cgi.K_dd_23_varindex,
+						      cgi.K_dd_33_varindex,
+	  cgi.psi_varindex,
+	  };
 const CCTK_INT input_array_type_codes[]
 	= {
 	  // g_ij
@@ -528,22 +528,53 @@ if (par_table_setup && (psi_flag == par_table_psi_flag))
 //
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling interpolator (%d points)",
-		   N_interp_points);
-status = CCTK_InterpLocalUniform(N_GRID_DIMS,
-				 gi.operator_handle, gi.param_table_handle,
-				 cgi.coord_origin, cgi.coord_delta,
-				 N_interp_points,
-				    interp_coords_type_code,
-				    interp_coords,
-				 N_input_arrays_use,
-				    cgi.gridfn_dims,
-				    input_array_type_codes,
-				    input_arrays,
-				 N_output_arrays_use,
-				    output_array_type_codes,
-				    output_arrays);
-
+		   "         calling %s interpolator (%d points)",
+		   global_local_str, 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.param_table_handle,
+				       gi.operator_handle,
+				       cgi.coord_system_handle,
+				       N_interp_points,
+					  interp_coords_type_code,
+					  interp_coords,
+				       N_input_arrays_use,
+					  input_array_variable_indices,
+				       N_output_arrays_use,
+					  output_array_type_codes,
+					  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,
+					 cgi.coord_origin, cgi.coord_delta,
+					 N_interp_points,
+					    interp_coords_type_code,
+					    interp_coords,
+					 N_input_arrays_use,
+					    cgi.gridfn_dims,
+					    input_array_type_codes,
+					    input_arrays,
+					 N_output_arrays_use,
+					    output_array_type_codes,
+					    output_arrays);
+	break;
+default:
+	error_exit(ERROR_EXIT,
+"***** interpolate_geometry(): bad gi.geometry_method=(int)%d!",
+		   int(gi.geometry_method));			/*NOTREACHED*/
+	}
 
 //
 // ***** handle any interpolation errors *****
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index c9cb8bc..40103ef 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -4,12 +4,11 @@
 // <<<prototypes for functions local to this file>>>
 //
 // expansion_Jacobian - top-level driver to compute the Jacobian
-// decode_geometry_method - decode the geometry_method parameter
 ///
 /// expansion_Jacobian_NP - compute the Jacobian by numerical perturbation
-///
-/// expansion_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
+/// expansion_Jacobian_partial_SD - compute partial-deriv terms: symbolic diff
 /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
+/// expansion_Jacobian_dr_FD - sum d/dr terms (compute via FD) into Jacobian
 ///
 
 #include <stdio.h>
@@ -134,29 +133,6 @@ return true;						// *** NORMAL RETURN ***
 }
 
 //******************************************************************************
-
-//
-// This function decodes the  Jacobian_method  parameter (string) into
-// an internal enum for future use.
-//
-enum Jacobian_method
-  decode_Jacobian_method(const char Jacobian_method_string[])
-{
-if	(STRING_EQUAL(Jacobian_method_string,
-		      "numerical perturbation"))
-   then return Jacobian_method__numerical_perturb;
-else if (STRING_EQUAL(Jacobian_method_string,
-		      "symbolic differentiation with finite diff d/dr"))
-   then return Jacobian_method__symbolic_diff_with_FD_dr;
-else if (STRING_EQUAL(Jacobian_method_string,
-		      "symbolic differentiation"))
-   then return Jacobian_method__symbolic_diff;
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"decode_Jacobian_method(): unknown Jacobian_method_string=\"%s\"!",
-		   Jacobian_method_string);			/*NOTREACHED*/
-}
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -260,8 +236,6 @@ return true;						// *** NORMAL RETURN ***
 	  }
 
 //******************************************************************************
-//******************************************************************************
-//******************************************************************************
 
 //
 // This function computes the partial derivative terms in the Jacobian
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 8ba54f8..b5e06f7 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -15,7 +15,8 @@ enum	{ N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 };
 //
 enum	geometry_method
 	{
-	geometry__interpolate_from_Cactus_grid,
+	geometry__global_interp_from_Cactus_grid,
+	geometry__local_interp_from_Cactus_grid,
 	geometry__Schwarzschild_EF // no comma
 	};
 
@@ -34,30 +35,27 @@ enum	Jacobian_method
 
 //
 // This structure holds all the information we need about the Cactus grid
-// and gridfns.  At present we use the  CCTK_InterpLocalUniform()  local
-// interpolator to access the Cactus gridfns.  Alas, much of the information
-// in this structure is specific to that API, and (alas) will need changing
-// when we eventually switch to the  CCTK_InterpGridArrays()  API.
+// and gridfns in order to interpolate the geometry information.
 //
 struct	cactus_grid_info
 	{
 	cGH *GH;			// --> Cactus grid hierarchy
 
-	// Cactus coordinate system
-	fp coord_origin[N_GRID_DIMS];	// (x,y,z) of grid posn (0,0,0)
-	fp coord_delta[N_GRID_DIMS];	// (x,y,z) grid spacing
-
-	// dimensions of gridfn data, viewed as a 3-D array
-	// n.b. storage ordering is Fortran,
-	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
-	CCTK_INT gridfn_dims[N_GRID_DIMS];
-
 	// this describes the semantics of the Cactus gij gridfns:
 	// true ==> the Cactus gij are conformal values as per
 	//          CactusEinstein/StaticConformal and the  psi  gridfn
 	// false ==> the Cactus gij are the physical metric
 	bool Cactus_conformal_metric;
 
+	//
+	// stuff for doing a global interpolation via CCTK_InterpGridArrays()
+	// i.e. geometry_info.geometry_method
+	//	== geometry__global_interp_from_Cactus_grid
+	//
+
+	// Cactus coordinate system
+	int coord_system_handle;
+
 	// Cactus variable indices
 	int g_dd_11_varindex, g_dd_12_varindex, g_dd_13_varindex,
 			      g_dd_22_varindex, g_dd_23_varindex,
@@ -67,6 +65,21 @@ struct	cactus_grid_info
 						K_dd_33_varindex;
 	int psi_varindex;		// unused if !Cactus_conformal_metric
 
+	//
+	// stuff for doing a local interpolation via CCTK_InterpLocalUniform()
+	// i.e. geometry_info.geometry_method
+	//	== geometry__local_interp_from_Cactus_grid
+	//
+
+	// Cactus coordinate system
+	fp coord_origin[N_GRID_DIMS];	// (x,y,z) of grid posn (0,0,0)
+	fp coord_delta[N_GRID_DIMS];	// (x,y,z) grid spacing
+
+	// dimensions of gridfn data, viewed as a 3-D array
+	// n.b. storage ordering is Fortran,
+	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
+	CCTK_INT gridfn_dims[N_GRID_DIMS];
+
 	// --> Cactus gridfn data for grid posn (0,0,0)
 	const fp* g_dd_11_data;
 	const fp* g_dd_12_data;
@@ -149,8 +162,6 @@ bool expansion(patch_system& ps,
 	       bool Jacobian_flag = false,
 	       bool print_msg_flag = false,
 	       jtutil::norm<fp>* H_norms_ptr = NULL);
-enum geometry_method
-  decode_geometry_method(const char geometry_method_string[]);
 
 // expansion_Jacobian.cc
 // ... returns true for successful computation, false for failure
@@ -160,10 +171,15 @@ bool expansion_Jacobian(patch_system& ps,
 			const struct geometry_info& gi,
 			const struct Jacobian_info& Jacobian_info,
 			bool print_msg_flag);
-enum Jacobian_method
-  decode_Jacobian_method(const char Jacobian_method_string[]);
 
 // Schwarzschild_EF.cc
 void Schwarzschild_EF_geometry(patch_system& ps,
 			       const struct geometry_info& gi,
 			       bool msg_flag);
+
+// misc.cc
+enum geometry_method
+  decode_geometry_method(const char geometry_method_string[]);
+bool geometry_method_is_interp(enum geometry_method geometry_method);
+enum Jacobian_method
+  decode_Jacobian_method(const char Jacobian_method_string[]);
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index e93414f..bffea61 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -4,7 +4,8 @@
 # Source files in this directory
 SRCS = expansion.cc \
        expansion_Jacobian.cc \
-       Schwarzschild_EF.cc
+       Schwarzschild_EF.cc \
+       misc.cc
 
 # Subdirectories containing source files
 SUBDIRS =