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

4 files changed, 108 insertions, 86 deletions

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 =