Fixed the names of keys for the table options entries which define the

boundary treatment behaviour for the local interpolator.
Also make a query call to the local interpolator to test whether it can deal
with the bounding box coordinates of each processor's local patch. This catches
errors in the multiprocessor case where the number of ghostzones isn't
sufficient to let the local interpolator do its job properly.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@46 1c20744c-e24a-42ec-9533-f5004cb800e5

1 files changed, 121 insertions, 92 deletions

diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c
index d9e9abe..3fc2ab0 100644
--- a/src/InterpGridArrays.c
+++ b/src/InterpGridArrays.c
@@ -185,7 +185,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
                                    void *const output_arrays[])
 {
   int i, retval;
-  CCTK_REAL *origin_local, *delta_local;
+  CCTK_REAL *origin_local, *delta;
   CCTK_INT *input_array_dims, *input_array_types, *input_array_time_levels;
   const void **input_arrays;
   const char *coord_system_name;
@@ -201,8 +201,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   char **output_arrays_local;
   type_desc_t *this, *type_desc;
   CCTK_REAL *range_min, *range_max;
-  CCTK_REAL *origin_global, *delta_global;
-  CCTK_REAL *interp_coords_proc, *coords;
+  CCTK_REAL *origin_global;
+  CCTK_REAL *interp_coords_proc, *coords, *bbox_local;
+  const CCTK_REAL **bbox_interp_coords;
 #endif
 
 
@@ -265,7 +266,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
 
   /* allocate some temporary arrays */
   origin_local = malloc (2 * N_dims * sizeof (CCTK_REAL));
-  delta_local  = origin_local + N_dims;
+  delta        = origin_local + N_dims;
   input_arrays = malloc (N_input_arrays * sizeof (void *));
   input_array_dims = malloc ((N_dims + 2*N_input_arrays) * sizeof (CCTK_INT));
   input_array_types = input_array_dims + N_dims;
@@ -298,9 +299,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
      copy the integer dimension array into an CCTK_INT array */
   for (i = 0; i < N_dims; i++)
   {
-    CCTK_CoordLocalRange (GH, &origin_local[i], &delta_local[i], i + 1,
-                          NULL, coord_system_name);
-    delta_local[i] = (delta_local[i] - origin_local[i]) / extras->lnsize[i];
+    CCTK_CoordLocalRange (GH, &origin_local[i], &delta[i], i + 1, NULL,
+                          coord_system_name);
+    delta[i] = (delta[i] - origin_local[i]) / extras->lnsize[i];
     input_array_dims[i] = extras->lnsize[i];
   }
 
@@ -354,32 +355,12 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     input_array_types[i] = CCTK_VarTypeI (input_array_indices[i]);
   }
 
-#if 0
-  /* Check whether the number of ghostzones is sufficient for the local
-     interpolator to handle points on the processor's local bounding box.
-     This is done as a query, simply by setting N_output_arrays = 0. */
-  if (retval == 0)
-  {
-    for (i = 0; i < N_dims; i++)
-    {
-      bbox_interp_coords[i] = origin_local;
-    retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
-                                      param_table_handle,
-                                      origin_local, delta_local, 2,
-                                      interp_coords_type, bbox_interp_coords,
-                                      N_input_arrays, input_array_dims,
-                                      input_array_types, input_arrays,
-                                      0, output_array_types,
-                                      output_arrays);
-  }
-#endif
-
   /* single-processor case directly calls the local interpolator */
   if (retval == 0 && pughGH->nprocs == 1)
   {
     retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
                                       param_table_handle,
-                                      origin_local, delta_local, N_points,
+                                      origin_local, delta, N_points,
                                       interp_coords_type, interp_coords,
                                       N_input_arrays, input_array_dims,
                                       input_array_types, input_arrays,
@@ -501,22 +482,21 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   {
     myGH->whichproc[point] = -1;
   }
-  
+
   /* initialize N_points_from to 0 for counting it up in the following loop */
   memset (myGH->N_points_from, 0, pughGH->nprocs * sizeof (CCTK_INT));
 
   /* allocate the ranges for my local coordinates */
-  range_min = malloc (4 * N_dims * sizeof (CCTK_REAL));
-  range_max = range_min + 1*N_dims;
+  range_min     = malloc (5 * N_dims * sizeof (CCTK_REAL));
+  range_max     = range_min + 1*N_dims;
   origin_global = range_min + 2*N_dims;
-  delta_global  = range_min + 3*N_dims;
+  bbox_local    = range_min + 3*N_dims;
 
-  /* get the global origin and delta of the coordinate system */
+  /* get the global origin of the coordinate system */
   for (i = 0; i < N_dims; i++)
   {
-    CCTK_CoordRange (GH, &origin_global[i], &delta_global[i], i+1, NULL,
+    CCTK_CoordRange (GH, &origin_global[i], &range_max[i], i+1, NULL,
                      coord_system_name);
-    delta_global[i] = (delta_global[i]-origin_global[i]) / (extras->nsize[i]-1);
   }
 
   /* locate the points to interpolate at */
@@ -526,15 +506,25 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     {
       /* compute the coordinate ranges */
       /* TODO: use bbox instead -- but the bboxes of other processors
-	 are not known */
+         are not known */
       int const has_lower = extras->lb[proc][i] == 0;
       int const has_upper = extras->ub[proc][i] == GH->cctk_gsh[i]-1;
-      range_min[i] = origin_global[i] + (extras->lb[proc][i] + (!has_lower) * (extras->nghostzones[i]-0.5) - FUDGE)*delta_global[i];
-      range_max[i] = origin_global[i] + (extras->ub[proc][i] - (!has_upper) * (extras->nghostzones[i]-0.5) + FUDGE)*delta_global[i];
+      range_min[i] = origin_global[i] + (extras->lb[proc][i] - FUDGE +
+                     (!has_lower) * (extras->nghostzones[i]-0.5)) * delta[i];
+      range_max[i] = origin_global[i] + (extras->ub[proc][i] + FUDGE -
+                     (!has_upper) * (extras->nghostzones[i]-0.5)) * delta[i];
+
 #ifdef PUGHINTERP_DEBUG
       printf("processor %d: range_min[%d]=%g range_max[%d]=%g\n",
              proc, i, (double) range_min[i], i, (double) range_max[i]);
 #endif
+
+      /* save this processor's bbox coordinates */
+      if (proc == pughGH->myproc)
+      {
+        bbox_local[2*i + 0] = range_min[i];
+        bbox_local[2*i + 1] = range_max[i];
+      }
     }
 
     /* and now which point will be processed by what processor */
@@ -563,29 +553,62 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
       }
     }
   }
-  /* don't need this anymore */
-  free (range_min);
+
+  /* Check whether the number of ghostzones is sufficient for the local
+     interpolator to handle points on the processor's local bounding box.
+     This is done as a query, simply by passing all arguments as for the
+     real call to the local interpolator, except for the output_arrays[]
+     argument which contains NULL pointers. */
+  bbox_interp_coords = malloc (N_dims * sizeof (void *));
+  for (i = 0; i < N_dims; i++)
+  {
+    bbox_interp_coords[i] = &bbox_local[2 * i];
+  }
+  output_arrays_local = N_output_arrays > 0 ?
+                        calloc (N_output_arrays, sizeof (void *)) : NULL;
+
+  retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
+                                    param_table_handle,
+                                    origin_local, delta, 2,
+                                    CCTK_VARIABLE_REAL,
+                                    (const void *const *) bbox_interp_coords,
+                                    N_input_arrays, input_array_dims,
+                                    input_array_types, input_arrays,
+                                    N_output_arrays, output_array_types,
+                                    (void *const *) output_arrays_local);
 
   /* make sure that all points could be mapped onto a processor */
-  for (point = j = 0; point < N_points; point++)
+  if (retval == 0)
   {
-    if (myGH->whichproc[point] < 0)
+    for (point = 0; point < N_points; point++)
     {
-      msg = malloc (80 + N_dims*20);
-      sprintf (msg, "Unable to locate point %d [%f",
-               point, (double) ((const CCTK_REAL *) interp_coords[0])[point]);
-      for (i = 1; i < N_dims; i++)
+      if (myGH->whichproc[point] < 0)
       {
-        sprintf (msg, "%s %f",
-                 msg, (double) ((const CCTK_REAL *) interp_coords[i])[point]);
+        msg = malloc (80 + N_dims*20);
+        sprintf (msg, "Unable to locate point %d [%f",
+                 point, (double) ((const CCTK_REAL *) interp_coords[0])[point]);
+        for (i = 1; i < N_dims; i++)
+        {
+          sprintf (msg, "%s %f",
+                   msg, (double) ((const CCTK_REAL *) interp_coords[i])[point]);
+        }
+        sprintf (msg, "%s]", msg);
+        CCTK_WARN (1, msg);
+        free (msg);
+        retval = UTIL_ERROR_BAD_INPUT;
       }
-      sprintf (msg, "%s]", msg);
-      CCTK_WARN (1, msg);
-      free (msg);
-      j = 1;                       /* mark as error */
     }
   }
-  if (j)
+
+  /* don't need this anymore */
+  if (output_arrays_local)
+  {
+    free (output_arrays_local);
+  }
+  free (bbox_interp_coords);
+  free (range_min);
+
+  if (retval)
   {
     if (myGH->whichproc)
     {
@@ -597,7 +620,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     free (input_array_dims);
     free (type_desc);
 
-    return (UTIL_ERROR_BAD_INPUT);
+    return (retval);
   }
 
   /* Now we want to resolve the N_points_from[]. Currently this is
@@ -631,7 +654,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
      This means I send interp_coords[i][point] to whichproc[point].
      I have N_points_from[proc] to send to each processor.
      */
-  
+
   /* This is backwards in the broadcast location; the number of points
      we are getting is how many everyone else is sending to us,
      eg, N_points_to, not how many we get back from everyone else,
@@ -683,7 +706,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     myGH->senddispl[proc] = myGH->senddispl[proc-1] + myGH->sendcnt[proc-1];
     myGH->recvdispl[proc] = myGH->recvdispl[proc-1] + myGH->recvcnt[proc-1];
   }
- 
+
   /* Great, and now exchange the coordinates and collect the ones
      that I have to process in *interp_coords_local[] */
   coords = NULL;
@@ -738,7 +761,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
        the results in the intermediate local output arrays */
     retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
                                       param_table_handle,
-                                      origin_local, delta_local, N_points_local,
+                                      origin_local, delta, N_points_local,
                                       interp_coords_type,
                                       (const void **) interp_coords_local,
                                       N_input_arrays, input_array_dims,
@@ -844,7 +867,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     }
 
     /* go back from processor-sorted data to input-ordered data.
-       The creation of the indices array above makes this not so bad. */ 
+       The creation of the indices array above makes this not so bad. */
     this->buf = this->recvbuf;
     for (proc = offset = 0; proc < pughGH->nprocs; proc++)
     {
@@ -967,7 +990,7 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
 {
   int i, retval;
   CCTK_INT *N_boundary_points_to_omit;
-  CCTK_REAL *off_centering_tolerance, *extrapolation_tolerance;
+  CCTK_REAL *boundary_off_centering_tolerance,*boundary_extrapolation_tolerance;
 
 
   retval = 0;
@@ -980,13 +1003,13 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
 
   /* allocate table option arrays and initialize them to their defaults */
   N_boundary_points_to_omit = malloc (2 * N_dims * sizeof (CCTK_INT));
-  off_centering_tolerance = malloc (2 * 2 * N_dims * sizeof (CCTK_REAL));
-  extrapolation_tolerance = off_centering_tolerance + 2 * N_dims;
+  boundary_off_centering_tolerance = malloc (2*2 * N_dims * sizeof (CCTK_REAL));
+  boundary_extrapolation_tolerance = boundary_off_centering_tolerance +2*N_dims;
   for (i = 0; i < 2 * N_dims; i++)
   {
     N_boundary_points_to_omit[i] = 0;
-    off_centering_tolerance[i] = 0.0;
-    extrapolation_tolerance[i] = 1e-10;
+    boundary_off_centering_tolerance[i] = 999.0;
+    boundary_extrapolation_tolerance[i] = 1e-10;
   }
 
   /* read the 'N_boundary_points_to_omit' options array from the table
@@ -996,54 +1019,57 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
                              "N_boundary_points_to_omit");
   if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
   {
-    CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' must be of "
-                  "CCTK_INT datatype");
+    CCTK_WARN (1, "Options array with key 'N_boundary_points_to_omit' "
+                  "must be of CCTK_INT datatype");
     retval = UTIL_ERROR_BAD_INPUT;
   }
   else if (i >= 0 && i < 2 * N_dims)
   {
-    CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' has fewer "
-                  "elements than interpolation bounding box");
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Options array with key 'N_boundary_points_to_omit' "
+                "must have exactly 2 * N_dims (= %d) elements", N_dims);
     retval = UTIL_ERROR_BAD_INPUT;
   }
 
-  /* read the 'off_centering_tolerance' options array from the table 
+  /* read the 'boundary_off_centering_tolerance' options array from the table
      and verify its datatype and element count */
   i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
-                              off_centering_tolerance,
-                              "off_centering_tolerance");
+                              boundary_off_centering_tolerance,
+                              "boundary_off_centering_tolerance");
   if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
   {
-    CCTK_WARN (1, "Array with key 'off_centering_tolerance' must be of "
-                  "CCTK_REAL datatype");
+    CCTK_WARN (1, "Options array with key 'boundary_off_centering_tolerance' "
+                  "must be of CCTK_REAL datatype");
     retval = UTIL_ERROR_BAD_INPUT;
   }
   else if (i >= 0 && i < 2 * N_dims)
   {
-    CCTK_WARN (1, "Array with key 'off_centering_tolerance' has fewer elements "
-                  "than interpolation bounding box");
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Options array with key 'boundary_off_centering_tolerance' "
+                "must have exactly 2 * N_dims (= %d) elements", N_dims);
     retval = UTIL_ERROR_BAD_INPUT;
   }
 
-  /* read the 'extrapolation_tolerance' options array from the table
+  /* read the 'boundary_extrapolation_tolerance' options array from the table
      and verify its datatype and element count */
   i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
-                              extrapolation_tolerance, 
-                              "extrapolation_tolerance");
+                              boundary_extrapolation_tolerance,
+                              "boundary_extrapolation_tolerance");
   if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
   {
-    CCTK_WARN (1, "Array with key 'extrapolation_tolerance' must be of "
-                  "CCTK_REAL datatype");
+    CCTK_WARN (1, "Options array with key 'boundary_extrapolation_tolerance' "
+                  "must be of CCTK_REAL datatype");
     retval = UTIL_ERROR_BAD_INPUT;
   }
-  else if (i >= 0 && i < 2 * N_dims)
+  else if (i >= 0 && i != 2 * N_dims)
   {
-    CCTK_WARN (1, "Array with key 'extrapolation_tolerance' has fewer elements "
-                  "than interpolation bounding box");
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Options array with key 'boundary_extrapolation_tolerance' "
+                "must have exactly 2 * N_dims (= %d) elements", N_dims);
     retval = UTIL_ERROR_BAD_INPUT;
   }
 
-  /* read the 'input_array_time_levels' options array from the table 
+  /* read the 'input_array_time_levels' options array from the table
      and verify its datatype and element count */
   i = Util_TableGetIntArray (param_table_handle, N_input_arrays,
                              input_array_time_levels,
@@ -1055,14 +1081,15 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
   }
   else if (i < 0)
   {
-    CCTK_WARN (1, "Array with key 'input_array_time_levels' must be of "
+    CCTK_WARN (1, "Options array with key 'input_array_time_levels' must be of "
                   "CCTK_INT datatype");
     retval = UTIL_ERROR_BAD_INPUT;
   }
-  else if (i < N_input_arrays)
+  else if (i != N_input_arrays)
   {
-    CCTK_WARN (1, "Array with key 'input_array_time_levels' has fewer elements "
-                  "than total number of input arrays");
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Options array with key 'input_array_time_levels' must have "
+                "exactly N_input_arrays (= %d) elements", N_input_arrays);
     retval = UTIL_ERROR_BAD_INPUT;
   }
 
@@ -1072,8 +1099,8 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
     if (! GH->cctk_bbox[i])
     {
       N_boundary_points_to_omit[i] = 0;
-      off_centering_tolerance[i] = 0.0;
-      extrapolation_tolerance[i] = 0.0;
+      boundary_off_centering_tolerance[i] = 0.0;
+      boundary_extrapolation_tolerance[i] = 0.0;
     }
   }
   if (! retval)
@@ -1082,14 +1109,16 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
                            N_boundary_points_to_omit,
                            "N_boundary_points_to_omit");
     Util_TableSetRealArray (param_table_handle, 2 * N_dims,
-                            off_centering_tolerance, "off_centering_tolerance");
+                            boundary_off_centering_tolerance,
+                            "boundary_off_centering_tolerance");
     Util_TableSetRealArray (param_table_handle, 2 * N_dims,
-                            extrapolation_tolerance, "extrapolation_tolerance");
+                            boundary_extrapolation_tolerance,
+                            "boundary_extrapolation_tolerance");
   }
 
   /* free allocated resources */
   free (N_boundary_points_to_omit);
-  free (off_centering_tolerance);
+  free (boundary_off_centering_tolerance);
 
   return (retval);
 }