2 bug fixes

examples/grid/in.sph

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

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) {
@@ -782,32 +785,54 @@ 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;
+            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;
+            continue;
+          }
+        } else {
+          skip[i] = 1;
+          continue;
+        }
       }
 
       skip[i] = 0;

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);
@@ -417,8 +427,8 @@ void Grid2d::initialize()
   recv_remap = nullptr;
 
   // store info about Comm decomposition needed for remap operation
-  // two Grid instances will exist for duration of remap
+  //   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;
@@ -593,9 +603,6 @@ void Grid2d::extract_comm_info()
     else if (rcbone.dim == 1) rcbone.cut = inylo;
     MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR,
                   rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm);
-    printf("G2D ExtractCommInfo: 0 dim/cut: %d %d, 1 dim/cut: %d %d\n",
-           rcbinfo[0].dim,rcbinfo[0].cut,
-           rcbinfo[1].dim,rcbinfo[1].cut);
   }
 }
 
@@ -617,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);
 }
 
@@ -1041,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();
 }
 
@@ -1087,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);
@@ -1192,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);
@@ -1652,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
 
@@ -1803,7 +1810,6 @@ void Grid2d::box_drop_grid(int *box, int proclower, int procupper,
 
   if (proclower == procupper) {
     plist[np++] = proclower;
-    printf("  proc %d plist np %d plist %d\n",comm->me,np,plist[np-1]);
     return;
   }
 
@@ -1815,13 +1821,9 @@ void Grid2d::box_drop_grid(int *box, int proclower, int procupper,
   // dim = 0,1,2 dimension of cut
 
   int procmid = proclower + (procupper - proclower) / 2 + 1;
-
-  printf("G2D BoxDropGrid: proc %d proc lum: %d %d %d np %d\n",comm->me,proclower,procupper,procmid,np);
-  //printf("  dim %d cut %d\n",rcbinfo[procmid].dim,rcbinfo[procmid].cut);
-
   int dim = rcbinfo[procmid].dim;
   int cut = rcbinfo[procmid].cut;
-  
+
   if (box[2*dim] < cut) box_drop_grid(box,proclower,procmid-1,np,plist);
   if (box[2*dim+1] >= cut) box_drop_grid(box,procmid,procupper,np,plist);
 }

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