Include ghostzones in processor-local hyperslabs.

This allows to recover from chunked datafiles/checkpoints which was broken
before because of the missing ghostzones.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHSlab/trunk@88 10716dce-81a3-4424-a2c8-48026a0d3035

1 files changed, 45 insertions, 35 deletions

diff --git a/src/Mapping.c b/src/Mapping.c
index f9e3803..00badad 100644
--- a/src/Mapping.c
+++ b/src/Mapping.c
@@ -39,14 +39,6 @@ CCTK_FILEVERSION(CactusPUGH_PUGHSlab_Mapping_c)
 #endif
 #define MAX(x, y)             ((x) > (y) ? (x) : (y))
 
-/* shortcuts for the local/global start/endpoints on this processor */
-#define MY_LOCAL_SP(extras, istag, dim)    (extras->ownership[istag][0][dim])
-#define MY_LOCAL_EP(extras, istag, dim)    (extras->ownership[istag][1][dim])
-#define MY_GLOBAL_SP(extras, myproc, istag, dim)                              \
-        (extras->lb[myproc][dim] + MY_LOCAL_SP (extras, istag, dim))
-#define MY_GLOBAL_EP(extras, myproc, istag, dim)                              \
-        (extras->lb[myproc][dim] + MY_LOCAL_EP (extras, istag, dim))
-
 
 /********************************************************************
  ********************    Static Variables    ************************
@@ -137,9 +129,10 @@ static CCTK_INT DefineMapping (const cGH      *GH,
 {
   int vdim, dim, num_dirs;
   int retval;
-  int stagger_index;
   int myproc, nprocs;
   int i, j, k, npoints;
+  int my_local_startpoint, my_local_endpoint;
+  int my_global_startpoint, my_global_endpoint;
   hslab_mapping_t *mapping;
   const char *error_msg;
   const pGH *pughGH;              /* pointer to the current pGH */
@@ -441,17 +434,28 @@ static CCTK_INT DefineMapping (const cGH      *GH,
     /* compute this processor's global startpoint from the global ranges */
     for (vdim = 0; vdim < vinfo.dim; vdim++)
     {
-      stagger_index = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+      i = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+
+      my_local_startpoint = extras->ownership[i][0][vdim];
+      my_local_endpoint   = extras->ownership[i][1][vdim];
+      /* local hyperslabs include ghostzones, global ones don't */
+      if (! is_global_hyperslab)
+      {
+        my_local_startpoint -= extras->nghostzones[vdim];
+        my_local_endpoint   += extras->nghostzones[vdim];
+      }
+
+      my_global_startpoint = extras->lb[myproc][vdim] + my_local_startpoint;
+      my_global_endpoint   = extras->lb[myproc][vdim] + my_local_endpoint;
 
-      if (origin[vdim] < MY_GLOBAL_EP (extras, myproc, stagger_index, vdim))
+      if (origin[vdim] < my_global_endpoint)
       {
         mapping->global_startpoint[vdim] = origin[vdim];
-        if (origin[vdim] < MY_GLOBAL_SP (extras, myproc,stagger_index,vdim))
+        if (origin[vdim] < my_global_startpoint)
         {
-          npoints = (MY_GLOBAL_SP (extras, myproc, stagger_index, vdim)
-                     - origin[vdim]) / mapping->downsample[vdim];
-          if ((MY_GLOBAL_SP (extras, myproc, stagger_index, vdim)
-               - origin[vdim]) % mapping->downsample[vdim])
+          npoints = (my_global_startpoint - origin[vdim]) /
+                    mapping->downsample[vdim];
+          if ((my_global_startpoint - origin[vdim]) % mapping->downsample[vdim])
           {
             npoints++;
           }
@@ -468,11 +472,22 @@ static CCTK_INT DefineMapping (const cGH      *GH,
     mapping->totals = 1;
     for (vdim = dim = 0; vdim < vinfo.dim; vdim++)
     {
-      stagger_index = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+      i = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+
+      my_local_startpoint = extras->ownership[i][0][vdim];
+      my_local_endpoint   = extras->ownership[i][1][vdim];
+      /* local hyperslabs include ghostzones, global ones don't */
+      if (! is_global_hyperslab)
+      {
+        my_local_startpoint -= extras->nghostzones[vdim];
+        my_local_endpoint   += extras->nghostzones[vdim];
+      }
+
+      my_global_startpoint = extras->lb[myproc][vdim] + my_local_startpoint;
+      my_global_endpoint   = extras->lb[myproc][vdim] + my_local_endpoint;
 
       if (mapping->global_startpoint[vdim] >= 0 &&
-          mapping->global_startpoint[vdim] <  MY_GLOBAL_EP (extras, myproc,
-                                                            stagger_index,vdim))
+          mapping->global_startpoint[vdim] <  my_global_endpoint)
       {
         mapping->local_startpoint[vdim] = mapping->global_startpoint[vdim] -
                                           extras->lb[myproc][vdim];
@@ -482,13 +497,11 @@ static CCTK_INT DefineMapping (const cGH      *GH,
         mapping->local_startpoint[vdim] = -1;
       }
 
-      if (mapping->global_endpoint[vdim] > MY_GLOBAL_SP (extras, myproc,
-                                                         stagger_index, vdim))
+      if (mapping->global_endpoint[vdim] > my_global_startpoint)
       {
         mapping->local_endpoint[vdim] =
-          MIN (MY_LOCAL_EP (extras, stagger_index, vdim),
-                            mapping->global_endpoint[vdim] -
-                            extras->lb[myproc][vdim]);
+          MIN (my_local_endpoint,
+               mapping->global_endpoint[vdim] - extras->lb[myproc][vdim]);
       }
       else
       {
@@ -496,10 +509,10 @@ static CCTK_INT DefineMapping (const cGH      *GH,
       }
 
 #ifdef DEBUG
-      printf ("direction %d: global ranges [%d, %d), local ranges[%d, %d)\n",
-              vdim,
-              mapping->global_startpoint[vdim], mapping->global_endpoint[vdim],
-              mapping->local_startpoint[vdim], mapping->local_endpoint[vdim]);
+      fprintf (stderr, "direction %d: global ranges [%d, %d), local ranges[%d, "
+               "%d)\n", vdim,
+               mapping->global_startpoint[vdim], mapping->global_endpoint[vdim],
+               mapping->local_startpoint[vdim], mapping->local_endpoint[vdim]);
 #endif
 
       if (mapping->local_endpoint[vdim] < 0 ||
@@ -513,8 +526,8 @@ static CCTK_INT DefineMapping (const cGH      *GH,
       {
         /* compute the local size in each hyperslab dimension */
         mapping->local_hsize[dim] = (mapping->local_endpoint[vdim] -
-                                   mapping->local_startpoint[vdim]) /
-                                  mapping->downsample[vdim];
+                                     mapping->local_startpoint[vdim]) /
+                                    mapping->downsample[vdim];
         if ((mapping->local_endpoint[vdim] - mapping->local_startpoint[vdim]) %
             mapping->downsample[vdim])
         {
@@ -526,12 +539,9 @@ static CCTK_INT DefineMapping (const cGH      *GH,
     }
   } /* end of else branch for 'if (mapping->is_full_hyperslab)' */
 
-#if 0
-  fprintf (stderr, "proc %d: total number of local hyperslab data points: %d\n", CCTK_MyProc (GH), mapping->totals);
-CCTK_Barrier (GH);
-#endif
 #ifdef DEBUG
-  printf ("total number of hyperslab data points: %d\n", mapping->totals);
+  fprintf (stderr, "total number of hyperslab data points: %d\n",
+           mapping->totals);
 #endif
 
   /* calculate the offset of the local chunk into the global hyperslab */