Another fix to return the correct status in the return code of CCTK_InterpGridArrays() and in the "local_interpolator_status" option in the user-supplied

parameter table.


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

1 files changed, 70 insertions, 52 deletions

diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c
index 53a1d90..8abc4bc 100644
--- a/src/InterpGridArrays.c
+++ b/src/InterpGridArrays.c
@@ -198,8 +198,8 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   CCTK_REAL **interp_coords_local;
   int N_points_local, N_types;
   int j, point, proc, array, offset, type, table_handle;
-  CCTK_INT have_per_point_status, error_point_status;
-  CCTK_INT *per_point_status, *per_proc_status, *local_interp_status;
+  CCTK_INT error_point_status;
+  CCTK_INT *per_point_status, *proc_status, *overall_status;
   char *msg;
   char **output_arrays_local;
   type_desc_t *this, *type_desc;
@@ -634,6 +634,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     retval = CCTK_ERROR_INTERP_GHOST_SIZE_TOO_SMALL;
   }
 
+#if 0
   /* make sure that all points could be mapped onto a processor */
   if (retval == 0)
   {
@@ -659,6 +660,30 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
       }
     }
   }
+#else
+  if (! retval)
+  {
+    /* points which couldn't be mapped onto any processor are now simply
+       mapped onto the local one and also passed on to the local interpolator */
+    for (point = 0; point < N_points; point++)
+    {
+      if (myGH->whichproc[point] < 0)
+      {
+        myGH->whichproc[point] = pughGH->myproc;
+        myGH->N_points_from[pughGH->myproc]++;
+      }
+    }
+  }
+  else
+  {
+    /* flag an error condition on this processor's error code in the
+       N_points_from[] array which then gets communicated to all processors */
+    for (proc = 0; proc < pughGH->nprocs; proc++)
+    {
+      myGH->N_points_from[proc] = retval;
+    }
+  }
+#endif
 
   /* don't need this anymore */
   if (output_arrays_local)
@@ -669,16 +694,6 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   free (bbox_interp_coords);
   free (range_min);
 
-  /* flag an error condition on this processor's error code in the
-     N_points_from[] array which then gets communicated to all processors */
-  if (retval)
-  {
-    for (proc = 0; proc < pughGH->nprocs; proc++)
-    {
-      myGH->N_points_from[proc] = retval;
-    }
-  }
-
   /* Now we want to resolve the N_points_from[]. Currently this is
      the form of ( in 2 proc mode )
      P1:  Num from P1  NFP2
@@ -835,18 +850,20 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     free (coords);
   }
 
-  /* prepare the per_proc_status[_p_] status array which will contain the
-     per-point status codes from CCTK_InterpLocalUniform() reduced over all
-     points requested by processor _p_
+  /* Prepare the proc_status[_p_] status array which in the first [nprocs]
+     elements will contain the minimum over all per-point status codes as
+     returned by CCTK_InterpLocalUniform() for all points requested by
+     processor _p_.
      Since the per_point_status is an optional return argument which may not
-     be supported by the local interpolator, we allocate the per_proc_status[]
-     array to be all zeros as a default status. */
-  have_per_point_status = 0;
-  per_proc_status = malloc ((2*pughGH->nprocs + N_points_local)
-                            * sizeof (CCTK_INT));
-  memset (per_proc_status, 0, pughGH->nprocs * sizeof (CCTK_INT));
-  local_interp_status = per_proc_status + 1*pughGH->nprocs;
-  per_point_status    = per_proc_status + 2*pughGH->nprocs;
+     be supported by the local interpolator, we allocate the proc_status[]
+     array to be all zeros as a default status.
+     The last element of the proc_status[] array will keep the return code
+     of CCTK_InterpLocalUniform() on processor _p_. */
+  proc_status = malloc ((2*(pughGH->nprocs + 1) + N_points_local)
+                        * sizeof (CCTK_INT));
+  memset (proc_status, 0, (pughGH->nprocs + 1) * sizeof (CCTK_INT));
+  overall_status   = proc_status + 1*(pughGH->nprocs + 1);
+  per_point_status = proc_status + 2*(pughGH->nprocs + 1);
 
   /* allocate intermediate output arrays for local interpolation */
   output_arrays_local = calloc (N_output_arrays, sizeof (void *));
@@ -876,33 +893,33 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
 
   /* now call the local interpolator for all local points and store
      the results in the intermediate local output arrays */
-  retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle, table_handle,
-                                    origin_local, delta, N_points_local,
-                                    interp_coords_type,
-                                    (const void **) interp_coords_local,
-                                    N_input_arrays, input_array_dims,
-                                    input_array_types, input_arrays,
-                                    N_output_arrays, output_array_types,
-                                    (void **) output_arrays_local);
+  proc_status[pughGH->nprocs] =
+    CCTK_InterpLocalUniform (N_dims, local_interp_handle, table_handle,
+                             origin_local, delta, N_points_local,
+                             interp_coords_type,
+                             (const void **) interp_coords_local,
+                             N_input_arrays, input_array_dims,
+                             input_array_types, input_arrays,
+                             N_output_arrays, output_array_types,
+                             (void **) output_arrays_local);
 
   /* Now check whether the local interpolator returned per-point
      status information.
      If so, and there were per-point errors returned by the local
      interpolator, then reduce the status values of all points from
      a requesting processor into a single value (ie. reduce
-     per_point_status[N_points_local] into per_proc_status[nprocs]).
-     The minimum overall status value will be taken. */
-  have_per_point_status = Util_TableGetInt (table_handle, &error_point_status,
-                                            "error_point_status") == 1;
-  if (have_per_point_status && error_point_status)
+     per_point_status[N_points_local] into proc_status[nprocs]).
+     Per definition the minimum overall status value will be taken. */
+  i = Util_TableGetInt (table_handle, &error_point_status,"error_point_status");
+  if (i == 1 && error_point_status)
   {
     for (proc = 0; proc < pughGH->nprocs; proc++)
     {
       for (point = 0; point < myGH->N_points_to[proc]; point++)
       {
-        if (per_proc_status[proc] > *per_point_status)
+        if (proc_status[proc] > *per_point_status)
         {
-          per_proc_status[proc] = *per_point_status;
+          proc_status[proc] = *per_point_status;
         }
         per_point_status++;
       }
@@ -920,17 +937,23 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     Util_TableDestroy (table_handle);
   }
 
-  /* If the local interpolator provided per-point status information
-     then do the global reduction on all per-processor status values now.
+  /* Now do the global reduction of all per-processor status values.
      This will then be taken as the return code to the caller. */
-  if (have_per_point_status)
+  CACTUS_MPI_ERROR (MPI_Allreduce (proc_status, overall_status,
+                                   pughGH->nprocs + 1, PUGH_MPI_INT, MPI_MIN,
+                                   pughGH->PUGH_COMM_WORLD));
+
+  /* if the local interpolator returned per-point status information then
+     store this processor's local interpolator status in the parameter table */
+  if (i == 1 && param_table_handle >= 0)
   {
-    CACTUS_MPI_ERROR (MPI_Allreduce (per_proc_status, local_interp_status,
-                                     pughGH->nprocs, PUGH_MPI_INT, MPI_MIN,
-                                     pughGH->PUGH_COMM_WORLD));
-    retval = local_interp_status[pughGH->myproc];
+    Util_TableSetInt (param_table_handle, overall_status[pughGH->myproc],
+                      "local_interpolator_status");
   }
 
+  /* set the return code to the overall processors' local interpolator status */
+  retval = overall_status[pughGH->nprocs];
+
   /* clean up some intermediate arrays */
   if (N_points_local > 0)
   {
@@ -939,7 +962,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
       free (interp_coords_local[i]);
     }
   }
-  free (per_proc_status);
+  free (proc_status);
   free (interp_coords_local);
   free (origin_local);
   free (input_arrays);
@@ -1099,12 +1122,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   free (type_desc);
 #endif /* MPI */
 
-  /* store this processor's interpolator status in the parameter table */
-  if (param_table_handle >= 0)
-  {
-    Util_TableSetInt (param_table_handle, retval, "local_interpolator_status");
-  }
-  return (0);
+  return (retval);
 }