Merge pull request #4289 from lammps/grid-debug

Fix 2 bugs with distributed grids

examples/grid/in.sph

@@ -109,4 +109,3 @@ variable           skin equal 0.3*${h}
 neighbor           ${skin} bin
 
 run                6000
-

src/compute_property_grid.cpp

@@ -255,8 +255,8 @@ void ComputePropertyGrid::allocate_grid()
 
   } else {
     grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
-    grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
+    grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in, nxlo_out, nxhi_out,
-                       nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
+                       nylo_out, nyhi_out, nzlo_out, nzhi_out);
     if (nvalues == 1)
       memory->create3d_offset(vec3d, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
                               "property/grid:vec3d");
@@ -342,13 +342,11 @@ void ComputePropertyGrid::pack_proc(int n)
     if (nvalues == 1) {
       for (int iz = nzlo_in; iz <= nzhi_in; iz++)
         for (int iy = nylo_in; iy <= nyhi_in; iy++)
-          for (int ix = nxlo_in; ix <= nxhi_in; ix++)
+          for (int ix = nxlo_in; ix <= nxhi_in; ix++) vec3d[iz][iy][ix] = me;
-            vec3d[iz][iy][ix] = me;
     } else {
       for (int iz = nzlo_in; iz <= nzhi_in; iz++)
         for (int iy = nylo_in; iy <= nyhi_in; iy++)
-          for (int ix = nxlo_in; ix <= nxhi_in; ix++)
+          for (int ix = nxlo_in; ix <= nxhi_in; ix++) array3d[iz][iy][ix][n] = me;
-            array3d[iz][iy][ix][n] = me;
     }
   }
 }
@@ -462,11 +460,15 @@ template <int POS, int MODE, int IDIM> void ComputePropertyGrid::pack_coords(int
 
       if (nvalues == 1) {
         for (int iy = nylo_in; iy <= nyhi_in; iy++) {
-          if (POS == LOW) lamda[1] = iy * dy;
+          if (POS == LOW)
-          else lamda[1] = (iy + 0.5) * dy;
+            lamda[1] = iy * dy;
+          else
+            lamda[1] = (iy + 0.5) * dy;
           for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
-            if (POS == LOW) lamda[0] = ix * dx;
+            if (POS == LOW)
-            else lamda[0] = (ix + 0.5) * dx;
+              lamda[0] = ix * dx;
+            else
+              lamda[0] = (ix + 0.5) * dx;
             domain->lamda2x(lamda, xone);
             if (IDIM == 0) vec2d[iy][ix] = xone[0];
             if (IDIM == 1) vec2d[iy][ix] = xone[1];
@@ -475,11 +477,15 @@ template <int POS, int MODE, int IDIM> void ComputePropertyGrid::pack_coords(int
 
       } else {
         for (int iy = nylo_in; iy <= nyhi_in; iy++) {
-          if (POS == LOW) lamda[1] = iy * dy;
+          if (POS == LOW)
-          else lamda[1] = (iy + 0.5) * dy;
+            lamda[1] = iy * dy;
+          else
+            lamda[1] = (iy + 0.5) * dy;
           for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
-            if (POS == LOW) lamda[0] = ix * dx;
+            if (POS == LOW)
-            else lamda[0] = (ix + 0.5) * dx;
+              lamda[0] = ix * dx;
+            else
+              lamda[0] = (ix + 0.5) * dx;
             domain->lamda2x(lamda, xone);
             if (IDIM == 0) array2d[iy][ix][n] = xone[0];
             if (IDIM == 1) array2d[iy][ix][n] = xone[1];
@@ -552,14 +558,20 @@ template <int POS, int MODE, int IDIM> void ComputePropertyGrid::pack_coords(int
 
       if (nvalues == 1) {
         for (int iz = nzlo_in; iz <= nzhi_in; iz++) {
-          if (POS == LOW) lamda[2] = iz * dz;
+          if (POS == LOW)
-          else lamda[2] = (iz + 0.5) * dz;
+            lamda[2] = iz * dz;
+          else
+            lamda[2] = (iz + 0.5) * dz;
           for (int iy = nylo_in; iy <= nyhi_in; iy++) {
-            if (POS == LOW) lamda[1] = iy * dy;
+            if (POS == LOW)
-            else lamda[1] = (iy + 0.5) * dy;
+              lamda[1] = iy * dy;
+            else
+              lamda[1] = (iy + 0.5) * dy;
             for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
-              if (POS == LOW) lamda[0] = ix * dx;
+              if (POS == LOW)
-              else lamda[0] = (ix + 0.5) * dx;
+                lamda[0] = ix * dx;
+              else
+                lamda[0] = (ix + 0.5) * dx;
               domain->lamda2x(lamda, xone);
               if (IDIM == 0) vec3d[iz][iy][ix] = xone[0];
               if (IDIM == 1) vec3d[iz][iy][ix] = xone[1];
@@ -570,14 +582,20 @@ template <int POS, int MODE, int IDIM> void ComputePropertyGrid::pack_coords(int
 
       } else {
         for (int iz = nzlo_in; iz <= nzhi_in; iz++) {
-          if (POS == LOW) lamda[2] = iz * dz;
+          if (POS == LOW)
-          else lamda[2] = (iz + 0.5) * dz;
+            lamda[2] = iz * dz;
+          else
+            lamda[2] = (iz + 0.5) * dz;
           for (int iy = nylo_in; iy <= nyhi_in; iy++) {
-            if (POS == LOW) lamda[1] = iy * dy;
+            if (POS == LOW)
-            else lamda[1] = (iy + 0.5) * dy;
+              lamda[1] = iy * dy;
+            else
+              lamda[1] = (iy + 0.5) * dy;
             for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
-              if (POS == LOW) lamda[0] = ix * dx;
+              if (POS == LOW)
-              else lamda[0] = (ix + 0.5) * dx;
+                lamda[0] = ix * dx;
+              else
+                lamda[0] = (ix + 0.5) * dx;
               domain->lamda2x(lamda, xone);
               if (IDIM == 0) array3d[iz][iy][ix][n] = xone[0];
               if (IDIM == 1) array3d[iz][iy][ix][n] = xone[1];

src/fix_ave_grid.cpp

@@ -684,11 +684,13 @@ void FixAveGrid::atom2grid()
   double ***vec3d = grid_sample->vec3d;
   double ****array3d = grid_sample->array3d;
 
+  // scan my owned atoms before tallying to bins
   // bin[i][dim] = indices of bin each atom is in
+  // count[][] or count[][][] = count of atoms contributing to each bin
+  // error check if any atom is out of bounds for my local+ghost grid cells
   // skip atom if group mask does not match
-  // check if any atom is out of bounds for my local grid
+  // skip atom if out of bounds for a nonperiodic dim and discardflag = DISCARD
-  // for nonperiodic dim, remap atom to first/last bin if out of bounds
+  // if out of bounds for a nonperiodic dim and discardflag = KEEP, remap atom to first/last bin
-  // count atoms contributing to each bin
 
   double *boxlo,*prd;
   int *periodicity = domain->periodicity;
@@ -719,6 +721,7 @@ void FixAveGrid::atom2grid()
 
   if (triclinic) domain->x2lamda(nlocal);
 
+  // set to 1
   int flag = 0;
 
   if (dimension == 2) {
@@ -732,23 +735,37 @@ void FixAveGrid::atom2grid()
       iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + OFFSET) - OFFSET;
 
       if (ix < nxlo_out || ix > nxhi_out) {
-        if (periodicity[0]) flag = 1;
+        if (periodicity[0]) {
-        else if (discardflag == KEEP) {
+          flag = 1;
+          continue;
+        } else if (discardflag == KEEP) {
           if (ix < nxlo_out && nxlo_out == 0) ix = 0;
           else if (ix > nxhi_out && nxhi_out == nxgrid-1) ix = nxgrid-1;
-          else flag = 1;
+          else {
-        } else skip[i] = 1;
+            flag = 1;
             continue;
           }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
+      }
       if (iy < nylo_out || iy > nyhi_out) {
-        if (periodicity[1]) flag = 1;
+        if (periodicity[1]) {
-        else if (discardflag == KEEP) {
+          flag = 1;
+          continue;
+        } else if (discardflag == KEEP) {
           if (iy < nylo_out && nylo_out == 0) iy = 0;
           else if (iy > nyhi_out && nyhi_out == nygrid-1) iy = nygrid-1;
-          else flag = 1;
+          else {
-        } else skip[i] = 1;
+            flag = 1;
             continue;
           }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
+      }
 
       skip[i] = 0;
       count2d[iy][ix] += 1.0;
@@ -768,33 +785,55 @@ void FixAveGrid::atom2grid()
       iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + OFFSET) - OFFSET;
 
       if (ix < nxlo_out || ix > nxhi_out) {
-        if (periodicity[0]) flag = 1;
+        if (periodicity[0]) {
-        else if (discardflag == KEEP) {
+          flag = 1;
+          continue;
+        } else if (discardflag == KEEP) {
           if (ix < nxlo_out && nxlo_out == 0) ix = 0;
           else if (ix > nxhi_out && nxhi_out == nxgrid-1) ix = nxgrid-1;
-          else flag = 1;
+          else {
-        } else skip[i] = 1;
+            flag = 1;
+            continue;
+          }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
       }
 
       if (iy < nylo_out || iy > nyhi_out) {
-        if (periodicity[1]) flag = 1;
+        if (periodicity[1]) {
-        else if (discardflag == KEEP) {
+          flag = 1;
+          continue;
+        } else if (discardflag == KEEP) {
           if (iy < nylo_out && nylo_out == 0) iy = 0;
           else if (iy > nyhi_out && nyhi_out == nygrid-1) iy = nygrid-1;
-          else flag = 1;
+          else {
-        } else skip[i] = 1;
+            flag = 1;
             continue;
           }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
+      }
 
       if (iz < nzlo_out || iz > nzhi_out) {
-        if (periodicity[2]) flag = 1;
+        if (periodicity[2]) {
-        else if (discardflag == KEEP) {
+          flag = 1;
+          continue;
+        } else if (discardflag == KEEP) {
           if (iz < nzlo_out && nzlo_out == 0) iz = 0;
           else if (iz > nzhi_out && nzhi_out == nzgrid-1) iz = nzgrid-1;
-          else flag = 1;
+          else {
-        } else skip[i] = 1;
+            flag = 1;
             continue;
           }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
+      }
 
       skip[i] = 0;
       count3d[iz][iy][ix] += 1.0;

src/fix_ave_grid.h

@@ -48,7 +48,7 @@ class FixAveGrid : public Fix {
   double memory_usage() override;
 
  private:
-  int nxgrid,nygrid,nzgrid;
+  int nxgrid, nygrid, nzgrid;
   int nvalues;
   int nrepeat, irepeat;
   bigint nvalid, nvalid_last;
@@ -57,14 +57,14 @@ class FixAveGrid : public Fix {
 
   double maxdist;
   int running_count;
-  int window_count,window_oldest,window_newest;
+  int window_count, window_oldest, window_newest;
 
   int biasflag;
   char *id_bias;
   class Compute *tbias;    // ptr to additional bias compute
-  double adof,cdof;
+  double adof, cdof;
 
-  int dimension,triclinic;
+  int dimension, triclinic;
 
   int *which, *argindex;
   char **ids;
@@ -75,26 +75,26 @@ class FixAveGrid : public Fix {
   int ngrid_buf1, ngrid_buf2;
   double *grid_buf1, *grid_buf2;
 
-  int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
+  int nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in;
-  int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
+  int nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out;
   int ngridout;
 
   struct GridData {
-    double **vec2d,***vec3d;
+    double **vec2d, ***vec3d;
-    double ***array2d,****array3d;
+    double ***array2d, ****array3d;
-    double **count2d,***count3d;
+    double **count2d, ***count3d;
   };
 
   GridData *grid_output;
-  GridData *grid_sample,*grid_nfreq,*grid_running;
+  GridData *grid_sample, *grid_nfreq, *grid_running;
   GridData **grid_window;
 
   // old grid data for remap operation
 
   class Grid2d *grid2d_previous;
   class Grid3d *grid3d_previous;
-  int nxlo_out_previous,nylo_out_previous,nzlo_out_previous;
+  int nxlo_out_previous, nylo_out_previous, nzlo_out_previous;
-  GridData *grid_sample_previous,*grid_nfreq_previous,*grid_running_previous;
+  GridData *grid_sample_previous, *grid_nfreq_previous, *grid_running_previous;
   GridData **grid_window_previous;
 
   int **bin;

src/grid2d.cpp

@@ -71,6 +71,11 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) :
   shift_atom_lo = shift_atom_hi = 0.0;
   yextra = 0;
   yfactor = 1.0;
+
+  // layout_grid = how this grid instance is distributed across procs
+  // depends on comm->layout at time this Grid2d instance is created
+
+  layout_grid = comm->layout;
 }
 
 /* ----------------------------------------------------------------------
@@ -112,6 +117,11 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int ixlo, int ixhi
   outylo = oylo;
   outyhi = oyhi;
 
+  // layout_grid = how this grid instance is distributed across procs
+  // depends on comm->layout at time this Grid2d instance is created
+
+  layout_grid = comm->layout;
+
   // additional intialization
   // other constructor invokes this from setup_grid()
 
@@ -338,7 +348,7 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, int &oxlo, i
 
   double fraclo, frachi;
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     fraclo = comm->xsplit[comm->myloc[0]];
     frachi = comm->xsplit[comm->myloc[0] + 1];
     partition_grid(nx, fraclo, frachi, shift_grid, 0, inxlo, inxhi);
@@ -418,7 +428,7 @@ void Grid2d::initialize()
 
   // store info about Comm decomposition needed for remap operation
   //   two Grid instances will exist for duration of remap
-  // each must know Comm decomp at time Grid instance was created
+  //   new Grid instance will then replace old one
 
   extract_comm_info();
 }
@@ -523,7 +533,7 @@ void Grid2d::ghost_grid()
   // also ensure no other procs use ghost cells beyond +y limit
 
   if (yextra) {
-    if (comm->layout != Comm::LAYOUT_TILED) {
+    if (layout_grid != Comm::LAYOUT_TILED) {
       if (comm->myloc[1] == comm->procgrid[1]-1) inyhi = outyhi = ny - 1;
     } else {
       if (comm->mysplit[1][1] == 1.0) inyhi = outyhi = ny - 1;
@@ -560,7 +570,7 @@ void Grid2d::extract_comm_info()
   // xyz split = copy of 1d vectors in Comm
   // grid2proc = copy of 3d array in Comm
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     procxlo = comm->procneigh[0][0];
     procxhi = comm->procneigh[0][1];
     procylo = comm->procneigh[1][0];
@@ -584,7 +594,7 @@ void Grid2d::extract_comm_info()
   // RCBinfo.cut = this proc's inlo in that dim
   // Allgather creates the tree of dims and cuts
 
-  if (comm->layout == Comm::LAYOUT_TILED) {
+  if (layout_grid == Comm::LAYOUT_TILED) {
     rcbinfo = (RCBinfo *)
       memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo");
     RCBinfo rcbone;
@@ -614,7 +624,7 @@ void Grid2d::extract_comm_info()
 
 void Grid2d::setup_comm(int &nbuf1, int &nbuf2)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2);
+  if (layout_grid != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2);
   else setup_comm_tiled(nbuf1,nbuf2);
 }
 
@@ -1038,7 +1048,7 @@ void Grid2d::setup_comm_tiled(int &nbuf1, int &nbuf2)
 
 int Grid2d::ghost_adjacent()
 {
-  if (comm->layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick();
+  if (layout_grid != Comm::LAYOUT_TILED) return ghost_adjacent_brick();
   return ghost_adjacent_tiled();
 }
 
@@ -1084,7 +1094,7 @@ int Grid2d::ghost_adjacent_tiled()
 void Grid2d::forward_comm(int caller, void *ptr, int which, int nper, int nbyte,
                           void *buf1, void *buf2, MPI_Datatype datatype)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     if (caller == KSPACE)
       forward_comm_brick<KSpace>((KSpace *) ptr,which,nper,nbyte,
                                  buf1,buf2,datatype);
@@ -1189,7 +1199,7 @@ forward_comm_tiled(T *ptr, int which, int nper, int nbyte,
 void Grid2d::reverse_comm(int caller, void *ptr, int which, int nper, int nbyte,
                           void *buf1, void *buf2, MPI_Datatype datatype)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     if (caller == KSPACE)
       reverse_comm_brick<KSpace>((KSpace *) ptr,which,nper,nbyte,
                                  buf1,buf2,datatype);
@@ -1649,7 +1659,7 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
 
   // test obox against appropriate layout
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
 
     // find comm->procgrid indices in each dim for box bounds
 
@@ -1713,7 +1723,7 @@ void Grid2d::clean_overlap()
 /* ----------------------------------------------------------------------
    recursively split a box until it doesn't overlap any periodic boundaries
    box = 4 integers = (xlo,xhi,ylo,yhi)
-     each lo/hi value may extend beyonw 0 to N-1 into another periodic image
+     each lo/hi value may extend beyond 0 to N-1 into another periodic image
    pbc = flags in each dim of which periodic image the caller box was in
    when a box straddles a periodic bounadry, split it in two
    when a box does not straddle, drop it down RCB tree

src/grid2d.h

@@ -57,6 +57,9 @@ class Grid2d : protected Pointers {
   int me, nprocs;
   MPI_Comm gridcomm;    // communicator for this class
                         // usually world, but MSM calls with subset
+  int layout_grid;      // how this grid instance is distributed across procs
+                        //   uses enum options for comm->layout
+                        //   load balancing can create a new Grid with new layout_grid
 
   // inputs from caller via constructor
 

src/grid3d.cpp

@@ -74,6 +74,11 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) :
   shift_atom_lo = shift_atom_hi = 0.0;
   zextra = 0;
   zfactor = 1.0;
+
+  // layout_grid = how this grid instance is distributed across procs
+  // depends on comm->layout at time this Grid3d instance is created
+
+  layout_grid = comm->layout;
 }
 
 /* ----------------------------------------------------------------------
@@ -122,6 +127,11 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz,
   outzlo = ozlo;
   outzhi = ozhi;
 
+  // layout_grid = how this grid instance is distributed across procs
+  // depends on comm->layout at time this Grid3d instance is created
+
+  layout_grid = comm->layout;
+
   // additional intialization
   // other constructor invokes this from setup_grid()
 
@@ -375,7 +385,7 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
 
   double fraclo,frachi;
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     fraclo = comm->xsplit[comm->myloc[0]];
     frachi = comm->xsplit[comm->myloc[0]+1];
     partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi);
@@ -579,7 +589,7 @@ void Grid3d::ghost_grid()
   // also ensure no other procs use ghost cells beyond +z limit
 
   if (zextra) {
-    if (comm->layout != Comm::LAYOUT_TILED) {
+    if (layout_grid != Comm::LAYOUT_TILED) {
       if (comm->myloc[2] == comm->procgrid[2]-1) inzhi = outzhi = nz - 1;
     } else {
       if (comm->mysplit[2][1] == 1.0) inzhi = outzhi = nz - 1;
@@ -621,7 +631,7 @@ void Grid3d::extract_comm_info()
   // xyz split = copy of 1d vectors in Comm
   // grid2proc = copy of 3d array in Comm
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     procxlo = comm->procneigh[0][0];
     procxhi = comm->procneigh[0][1];
     procylo = comm->procneigh[1][0];
@@ -649,7 +659,7 @@ void Grid3d::extract_comm_info()
   // RCBinfo.cut = this proc's inlo in that dim
   // Allgather creates the tree of dims and cuts
 
-  if (comm->layout == Comm::LAYOUT_TILED) {
+  if (layout_grid == Comm::LAYOUT_TILED) {
     rcbinfo = (RCBinfo *)
       memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo");
     RCBinfo rcbone;
@@ -680,7 +690,7 @@ void Grid3d::extract_comm_info()
 
 void Grid3d::setup_comm(int &nbuf1, int &nbuf2)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2);
+  if (layout_grid != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2);
   else setup_comm_tiled(nbuf1,nbuf2);
 }
 
@@ -1207,7 +1217,7 @@ void Grid3d::setup_comm_tiled(int &nbuf1, int &nbuf2)
 
 int Grid3d::ghost_adjacent()
 {
-  if (comm->layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick();
+  if (layout_grid != Comm::LAYOUT_TILED) return ghost_adjacent_brick();
   return ghost_adjacent_tiled();
 }
 
@@ -1255,7 +1265,7 @@ int Grid3d::ghost_adjacent_tiled()
 void Grid3d::forward_comm(int caller, void *ptr, int which, int nper, int nbyte,
                             void *buf1, void *buf2, MPI_Datatype datatype)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     if (caller == KSPACE)
       forward_comm_brick<KSpace>((KSpace *) ptr,which,nper,nbyte,
                                  buf1,buf2,datatype);
@@ -1360,7 +1370,7 @@ forward_comm_tiled(T *ptr, int which, int nper, int nbyte,
 void Grid3d::reverse_comm(int caller, void *ptr, int which, int nper, int nbyte,
                             void *buf1, void *buf2, MPI_Datatype datatype)
 {
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
     if (caller == KSPACE)
       reverse_comm_brick<KSpace>((KSpace *) ptr,which,nper,nbyte,
                                 buf1,buf2,datatype);
@@ -1825,7 +1835,7 @@ int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
     return noverlap_list;
   }
 
-  if (comm->layout != Comm::LAYOUT_TILED) {
+  if (layout_grid != Comm::LAYOUT_TILED) {
 
     // find comm->procgrid indices in each dim for box bounds
 

src/grid3d.h

@@ -59,6 +59,9 @@ class Grid3d : protected Pointers {
   int me, nprocs;
   MPI_Comm gridcomm;    // communicator for this class
                         // usually world, but MSM calls with subset
+  int layout_grid;      // how this grid instance is distributed across procs
+                        //   uses enum options for comm->layout
+                        //   load balancing can create a new Grid with new layout_grid
 
   // input from caller