bug fixes to grid remapping

src/dump_grid.cpp

@@ -26,6 +26,9 @@
 #include "region.h"
 #include "update.h"
 
+// DEBUG
+#include "comm.h"
+
 #include <cstring>
 
 using namespace LAMMPS_NS;

src/fix_ave_grid.cpp

@@ -29,9 +29,6 @@
 #include "update.h"
 #include "variable.h"
 
-// DEBUG
-#include "comm.h"
-
 #include <cstring>
 
 using namespace LAMMPS_NS;
@@ -391,65 +388,24 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
-  // instantiate the Grid class and allocate per-grid memory
+  // instantiate Grid class and buffers
+  // allocate/zero per-grid data
 
-  if (modeatom) maxdist = 0.5 * neighbor->skin;
+  allocate_grid();
-  else if (modegrid) maxdist = 0.0;
 
-  if (dimension == 2) {
+  grid_sample = allocate_one_grid();
-    grid2d = new Grid2d(lmp, world, nxgrid, nygrid);
+  grid_nfreq = allocate_one_grid();
-    grid2d->set_distance(maxdist);
+  if (aveflag == RUNNING || aveflag == WINDOW) grid_nfreq = allocate_one_grid();
-    grid2d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in,
-                       nxlo_out, nxhi_out, nylo_out, nyhi_out);
-
-    // ngrid_buf12 converted to nvalues + count
-
-    grid2d->setup_comm(ngrid_buf1, ngrid_buf2);
-    ngrid_buf1 *= nvalues + 1;
-    ngrid_buf2 *= nvalues + 1;
-
-    memory->create(grid_buf1, ngrid_buf1, "ave/grid:grid_buf1");
-    memory->create(grid_buf2, ngrid_buf2, "ave/grid:grid_buf2");
-
-    ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1);
-
-  } else {
-    grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
-    grid3d->set_distance(maxdist);
-    grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
-                       nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
-
-    // ngrid_buf12 converted to nvalues + count
-
-    grid3d->setup_comm(ngrid_buf1, ngrid_buf2);
-    ngrid_buf1 *= nvalues + 1;
-    ngrid_buf2 *= nvalues + 1;
-
-    memory->create(grid_buf1, ngrid_buf1, "ave/grid:grid_buf1");
-    memory->create(grid_buf2, ngrid_buf2, "ave/grid:grid_buf2");
-
-    ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1) *
-      (nzhi_out - nzlo_out + 1);
-  }
-
-  // create data structs for per-grid data
-
-  grid_output = new GridData();
-  grid_sample = new GridData();
-  grid_nfreq = new GridData();
-  grid_running = new GridData();
   if (aveflag == WINDOW) {
     grid_window = new GridData*[nwindow];
     for (int i = 0; i < nwindow; i++)
-      grid_window[i] = new GridData();
+      grid_window[i] = allocate_one_grid();
-  } else grid_window = nullptr;
+  }
 
-  allocate_grid(grid_sample);
+  // output may occur via dump on timestep 0
-  allocate_grid(grid_nfreq);
+
-  if (aveflag == RUNNING || aveflag == WINDOW) allocate_grid(grid_running);
+  grid_output = new GridData();
-  if (aveflag == WINDOW)
+  output_grid(grid_nfreq);
-    for (int i = 0; i < nwindow; i++)
-      allocate_grid(grid_window[i]);
 
   // initialize running and window values
 
@@ -458,14 +414,6 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
   window_oldest = -1;
   window_newest = 0;
 
-  // zero grid_nfreq for output since dump may access it on timestep 0
-  // also one-time zero of grid_running for ave = RUNNING or WINDOW
-
-  zero_grid(grid_nfreq);
-  output_grid(grid_nfreq);
-
-  if (aveflag == RUNNING || aveflag == WINDOW) zero_grid(grid_running);
-
   // bin indices and skip flags for ATOM mode
   // vresult for per-atom variable evaluation
 
@@ -499,29 +447,28 @@ FixAveGrid::~FixAveGrid()
   delete[] value2grid;
   delete[] value2data;
 
-  delete grid2d;
+  // deallocate Grid class and buffers
-  delete grid3d;
+  
+  if (dimension == 2) delete grid2d;
+  else delete grid3d;
 
   memory->destroy(grid_buf1);
   memory->destroy(grid_buf2);
 
-  // deallocate all per-grid data
+  // deallocate per-grid data
 
-  deallocate_grid(grid_sample);
+  deallocate_one_grid(grid_sample,nxlo_out,nylo_out,nzlo_out);
-  deallocate_grid(grid_nfreq);
+  deallocate_one_grid(grid_nfreq,nxlo_out,nylo_out,nzlo_out);
-  if (aveflag == RUNNING || aveflag == WINDOW) deallocate_grid(grid_running);
+  if (aveflag == RUNNING || aveflag == WINDOW)
-  if (aveflag == WINDOW)
+    deallocate_one_grid(grid_running,nxlo_out,nylo_out,nzlo_out);
-    for (int i = 0; i < nwindow; i++) {
+  if (aveflag == WINDOW) {
-      deallocate_grid(grid_window[i]);
+    for (int i = 0; i < nwindow; i++)
-      delete grid_window[i];
+      deallocate_one_grid(grid_window[i],nxlo_out,nylo_out,nzlo_out);
-    }
+    delete [] grid_window;
+  }
 
   delete grid_output;
-  delete grid_sample;
+  
-  delete grid_nfreq;
-  delete grid_running;
-  delete [] grid_window;
-
   if (modeatom) {
     memory->destroy(bin);
     memory->destroy(skip);
@@ -719,12 +666,12 @@ void FixAveGrid::end_of_step()
   //   for norm = ALL, normalize sample grid by counts over all samples
   //   for norm = SAMPLE, normalize Nfreq grid by Nrepeat
   //   for norm = NONORM, normalize sample grid by Nrepeat, not by counts
-  //              this check is made inside normalize_grid()
+  //              this check is made inside normalize_atom()
   // for GRID mode:
   //   normalize sample grid by Nrepeat
 
   if (modeatom) {
-    if (normflag == ALL) {
+    if (normflag == ALL) { 
       normalize_atom(nrepeat,grid_sample);
       normalize_count(nrepeat,grid_sample);
       copy_grid(grid_sample,grid_nfreq);
@@ -1342,7 +1289,7 @@ void FixAveGrid::normalize_atom(int numsamples, GridData *grid)
             else if (which[0] == ArgInfo::DENSITY_MASS)
               norm = density_mass_norm;
             else if (which[0] == ArgInfo::TEMPERATURE)
-              norm = mvv2e /((repeat*cdof + adof*count) * boltz);
+              norm = mvv2e / ((repeat*cdof + adof*count) * boltz);
             else if (normflag == NONORM)
               norm = invrepeat;
             else
@@ -1363,7 +1310,7 @@ void FixAveGrid::normalize_atom(int numsamples, GridData *grid)
               else if (which[m] == ArgInfo::DENSITY_MASS)
                 norm = density_mass_norm;
               else if (which[m] == ArgInfo::TEMPERATURE)
-                norm = mvv2e /((repeat*cdof + adof*count) * boltz);
+                norm = mvv2e / ((repeat*cdof + adof*count) * boltz);
               else if (normflag == NONORM)
                 norm = invrepeat;
               else
@@ -1389,7 +1336,7 @@ void FixAveGrid::normalize_atom(int numsamples, GridData *grid)
               else if (which[0] == ArgInfo::DENSITY_MASS)
                 norm = density_mass_norm;
               else if (which[0] == ArgInfo::TEMPERATURE)
-                norm = mvv2e /((repeat*cdof + adof*count) * boltz);
+                norm = mvv2e / ((repeat*cdof + adof*count) * boltz);
               else if (normflag == NONORM)
                 norm = invrepeat;
               else
@@ -1411,7 +1358,7 @@ void FixAveGrid::normalize_atom(int numsamples, GridData *grid)
                 else if (which[m] == ArgInfo::DENSITY_MASS)
                 norm = density_mass_norm;
                 else if (which[m] == ArgInfo::TEMPERATURE)
-                  norm = mvv2e /((repeat*cdof + adof*count) * boltz);
+                  norm = mvv2e / ((repeat*cdof + adof*count) * boltz);
                 else if (normflag == NONORM)
                   norm = invrepeat;
                 else
@@ -1495,12 +1442,69 @@ void FixAveGrid::normalize_count(int numsamples, GridData *grid)
 }
 
 /* ----------------------------------------------------------------------
-   allocate a data grid
+   allocate instance of Grid2d or Grid3d
+------------------------------------------------------------------------- */
+
+void FixAveGrid::allocate_grid()
+{
+  if (modeatom) maxdist = 0.5 * neighbor->skin;
+  else if (modegrid) maxdist = 0.0;
+
+  if (dimension == 2) {
+    grid2d = new Grid2d(lmp, world, nxgrid, nygrid);
+    grid2d->set_distance(maxdist);
+    grid2d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in,
+                       nxlo_out, nxhi_out, nylo_out, nyhi_out);
+
+    // ngrid_buf12 converted to nvalues + count
+
+    grid2d->setup_comm(ngrid_buf1, ngrid_buf2);
+    ngrid_buf1 *= nvalues + 1;
+    ngrid_buf2 *= nvalues + 1;
+
+    grid_buf1 = grid_buf2 = nullptr;
+    if (ngrid_buf1) memory->create(grid_buf1, ngrid_buf1, "ave/grid:grid_buf1");
+    if (ngrid_buf2) memory->create(grid_buf2, ngrid_buf2, "ave/grid:grid_buf2");
+
+    ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1);
+
+  } else {
+    grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
+    grid3d->set_distance(maxdist);
+    grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
+                       nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
+
+    // ngrid_buf12 converted to nvalues + count
+
+    grid3d->setup_comm(ngrid_buf1, ngrid_buf2);
+    ngrid_buf1 *= nvalues + 1;
+    ngrid_buf2 *= nvalues + 1;
+
+    grid_buf1 = grid_buf2 = nullptr;
+    if (ngrid_buf1) memory->create(grid_buf1, ngrid_buf1, "ave/grid:grid_buf1");
+    if (ngrid_buf2) memory->create(grid_buf2, ngrid_buf2, "ave/grid:grid_buf2");
+
+    ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1) *
+      (nzhi_out - nzlo_out + 1);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   allocate a data grid and zero its values
    if ATOM mode, also allocate per-grid count
 ------------------------------------------------------------------------- */
 
-void FixAveGrid::allocate_grid(GridData *grid)
+FixAveGrid::GridData *FixAveGrid::allocate_one_grid()
 {
+  GridData *grid = new GridData();
+
+  grid->vec2d = nullptr;
+  grid->array2d = nullptr;
+  grid->count2d = nullptr;
+  grid->vec3d = nullptr;
+  grid->array3d = nullptr;
+  grid->count3d = nullptr;
+  
   if (dimension == 2) {
     if (nvalues == 1)
       memory->create2d_offset(grid->vec2d, nylo_out, nyhi_out,
@@ -1525,31 +1529,59 @@ void FixAveGrid::allocate_grid(GridData *grid)
       memory->create3d_offset(grid->count3d, nzlo_out, nzhi_out, nylo_out,
                               nyhi_out, nxlo_out, nxhi_out, "ave/grid:count3d");
   }
+
+  zero_grid(grid);
+  
+  return grid;
+}
+
+
+/* ----------------------------------------------------------------------
+   create clone of a data grid
+   allocate a new grid and copy only the pointers from the source grid
+   used by reset_grid() to keep old data values until remap is complete
+------------------------------------------------------------------------- */
+
+FixAveGrid::GridData *FixAveGrid::clone_one_grid(GridData *src)
+{
+  GridData *grid = new GridData();
+
+  grid->vec2d = src->vec2d;
+  grid->array2d = src->array2d;
+  grid->count2d = src->count2d;
+  grid->vec3d = src->vec3d;
+  grid->array3d = src->array3d;
+  grid->count3d = src->count3d;
+
+  return grid;
 }
 
 /* ----------------------------------------------------------------------
-   deallocate a data grid
+   deallocate a data grid and all its memory
    if ATOM mode, also deallocate per-grid count
 ------------------------------------------------------------------------- */
 
-void FixAveGrid::deallocate_grid(GridData *grid)
+void FixAveGrid::deallocate_one_grid(GridData *grid,
+                                     int xoffset, int yoffset, int zoffset)
 {
   if (dimension == 2) {
     if (nvalues == 1)
-      memory->destroy2d_offset(grid->vec2d,nylo_out,nxlo_out);
+      memory->destroy2d_offset(grid->vec2d,yoffset,xoffset);
     else
-      memory->destroy3d_offset_last(grid->array2d,nylo_out,nxlo_out);
+      memory->destroy3d_offset_last(grid->array2d,yoffset,xoffset);
     if (modeatom)
-      memory->destroy2d_offset(grid->count2d,nylo_out,nxlo_out);
+      memory->destroy2d_offset(grid->count2d,yoffset,xoffset);
 
   } else if (dimension == 3) {
     if (nvalues == 1)
-      memory->destroy3d_offset(grid->vec3d,nzlo_out,nylo_out,nxlo_out);
+      memory->destroy3d_offset(grid->vec3d,zoffset,yoffset,xoffset);
     else
-      memory->destroy4d_offset_last(grid->array3d,nzlo_out,nylo_out,nxlo_out);
+      memory->destroy4d_offset_last(grid->array3d,zoffset,yoffset,xoffset);
     if (modeatom)
-      memory->destroy3d_offset(grid->count3d,nzlo_out,nylo_out,nxlo_out);
+      memory->destroy3d_offset(grid->count3d,zoffset,yoffset,xoffset);
   }
+
+  delete grid;
 }
 
 /* ----------------------------------------------------------------------
@@ -1591,7 +1623,7 @@ void FixAveGrid::zero_grid(GridData *grid)
     if (modeatom)
       memset(&grid->count2d[nylo_out][nxlo_out],0,ngridout*sizeof(double));
 
-  } else if (dimension == 3) {
+  } else {
     if (nvalues == 1)
       memset(&grid->vec3d[nzlo_out][nylo_out][nxlo_out],0,
              ngridout*sizeof(double));
@@ -1817,8 +1849,9 @@ void FixAveGrid::output_grid(GridData *src)
 
 /* ----------------------------------------------------------------------
    pack ghost values into buf to send to another proc
-   nvalues per grid point + count
+   nvalues per grid point + count-
-------------------------------------------------------------------------- */
+   only invoked for ATOM mode
+------------------------------------------------------------------------ */
 
 void FixAveGrid::pack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *list)
 {
@@ -1837,7 +1870,7 @@ void FixAveGrid::pack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *lis
     if (nvalues == 1) data = &grid_sample->vec3d[nzlo_out][nylo_out][nxlo_out];
     else data = &grid_sample->array3d[nzlo_out][nylo_out][nxlo_out][0];
   }
-
+  
   if (nvalues == 1) {
     for (i = 0; i < nlist; i++) {
       buf[m++] = count[list[i]];
@@ -1856,6 +1889,7 @@ void FixAveGrid::pack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *lis
 /* ----------------------------------------------------------------------
    unpack another proc's ghost values from buf and add to own values
    nvalues per grid point + count
+   only invoked for ATOM mode
 ------------------------------------------------------------------------- */
 
 void FixAveGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *list)
@@ -1864,7 +1898,6 @@ void FixAveGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *l
 
   auto buf = (double *) vbuf;
   double *count,*data,*values;
-  m = 0;
 
   if (dimension == 2) {
     count = &grid_sample->count2d[nylo_out][nxlo_out];
@@ -1876,6 +1909,7 @@ void FixAveGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *l
     else data = &grid_sample->array3d[nzlo_out][nylo_out][nxlo_out][0];
   }
 
+  m = 0;
   if (nvalues == 1) {
     for (i = 0; i < nlist; i++) {
       count[list[i]] += buf[m++];
@@ -1893,29 +1927,118 @@ void FixAveGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *l
 
 /* ----------------------------------------------------------------------
    pack old grid values to buf to send to another proc
+   invoked for both GRID and ATOM mode
 ------------------------------------------------------------------------- */
 
 void FixAveGrid::pack_remap_grid(void *vbuf, int nlist, int *list)
 {
-  auto buf = (double *) vbuf;
+  int i,j,m,iwindow;
-  double *src;
-  //double *src =
-  //  &T_electron_previous[nzlo_out_previous][nylo_out_previous][nxlo_out_previous];
 
-  for (int i = 0; i < nlist; i++) buf[i] = src[list[i]];
+  auto buf = (double *) vbuf;
+
+  int running_flag = 0;
+  if (aveflag == RUNNING || aveflag == WINDOW) running_flag = 1;
+  int window_flag = 0;
+  if (aveflag == WINDOW) window_flag = 1;
+
+  m = 0;
+  for (i = 0; i < nlist; i++) { 
+    m += pack_one_grid(grid_sample_previous,list[i],&buf[m]);
+    m += pack_one_grid(grid_nfreq_previous,list[i],&buf[m]);
+    if (running_flag) m += pack_one_grid(grid_running_previous,list[i],&buf[m]);
+    if (window_flag)
+      for (iwindow = 0; iwindow < nwindow; iwindow++)
+        m += pack_one_grid(grid_window_previous[iwindow],list[i],&buf[m]);
+  }
 }
 
 /* ----------------------------------------------------------------------
-   unpack another proc's own values from buf and set own ghost values
+   unpack received owned values from buf into new grid
+   invoked for both GRID and ATOM mode
 ------------------------------------------------------------------------- */
 
 void FixAveGrid::unpack_remap_grid(void *vbuf, int nlist, int *list)
 {
-  auto buf = (double *) vbuf;
+  int i,j,m,iwindow;
-  double *dest;
-  //double *dest = &T_electron[nzlo_out][nylo_out][nxlo_out];
 
-  for (int i = 0; i < nlist; i++) dest[list[i]] = buf[i];
+  auto buf = (double *) vbuf;
+
+  int running_flag = 0;
+  if (aveflag == RUNNING || aveflag == WINDOW) running_flag = 1;
+  int window_flag = 0;
+  if (aveflag == WINDOW) window_flag = 1;
+
+  m = 0;
+  for (i = 0; i < nlist; i++) {
+    m += unpack_one_grid(&buf[m],grid_sample,list[i]);
+    m += unpack_one_grid(&buf[m],grid_nfreq,list[i]);
+    if (running_flag) m += unpack_one_grid(&buf[m],grid_running,list[i]);
+    if (window_flag)
+      for (iwindow = 0; iwindow < nwindow; iwindow++)
+        m += unpack_one_grid(&buf[m],grid_window[iwindow],list[i]);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   pack values for a single grid cell to buf
+   return number of values packed
+------------------------------------------------------------------------- */
+
+int FixAveGrid::pack_one_grid(GridData *grid, int index, double *buf)
+{
+  double *count,*data,*values;
+
+  if (dimension == 2) {
+    count = &grid->count2d[nylo_out_previous][nxlo_out_previous];
+    if (nvalues == 1) data = &grid->vec2d[nylo_out_previous][nxlo_out_previous];
+    else data = &grid->array2d[nylo_out_previous][nxlo_out_previous][0];
+  } else if (dimension == 3) {
+    count = &grid->count3d[nzlo_out_previous][nylo_out_previous][nxlo_out_previous];
+    if (nvalues == 1) data = &grid->vec3d[nzlo_out_previous][nylo_out_previous][nxlo_out_previous];
+    else data = &grid->array3d[nzlo_out_previous][nylo_out_previous][nxlo_out_previous][0];
+  }
+
+  int m = 0;
+  if (modeatom) buf[m++] = count[index];
+  if (nvalues == 1) buf[m++] = data[index];
+  else {
+    values = &data[nvalues*index];
+    for (int j = 0; j < nvalues; j++)
+      buf[m++] = values[j];
+  }
+
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values for a single grid cell from buf
+   return number of values unpacked
+------------------------------------------------------------------------- */
+
+int FixAveGrid::unpack_one_grid(double *buf, GridData *grid, int index)
+{
+  double *count,*data,*values;
+
+  if (dimension == 2) {
+    count = &grid->count2d[nylo_out][nxlo_out];
+    if (nvalues == 1) data = &grid->vec2d[nylo_out][nxlo_out];
+    else data = &grid->array2d[nylo_out][nxlo_out][0];
+  } else if (dimension == 3) {
+    count = &grid->count3d[nzlo_out][nylo_out][nxlo_out];
+    if (nvalues == 1) data = &grid->vec3d[nzlo_out][nylo_out][nxlo_out];
+    else data = &grid->array3d[nzlo_out][nylo_out][nxlo_out][0];
+  }
+
+  int m = 0;
+  if (modeatom) count[index] = buf[m++];
+  if (nvalues == 1) data[index] = buf[m++];
+  else {
+    values = &data[nvalues*index];
+    for (int j = 0; j < nvalues; j++)
+      values[j] = buf[m++];
+  }
+
+  return m;
 }
 
 /* ----------------------------------------------------------------------
@@ -1955,57 +2078,104 @@ void FixAveGrid::reset_grid()
     } else delete gridnew;
   }
 
-  // DEBUG
+  // for ATOM mode, perform ghost to owned grid comm for grid_sample
-  if (comm->me == 0) printf("Remapping grid on step %ld\n",update->ntimestep);
+  // necessary b/c remap will only communicate owned grid cell data
+  //   so can't lose ghost data not yet summed to owned cells
+  // nvalues+1 includes atom count
 
-  /*
+  if (modeatom) {
-  // delete grid data which doesn't need to persist from previous to new decomp
+    if (dimension == 2)
+      grid2d->reverse_comm(Grid2d::FIX,this,nvalues+1,sizeof(double),0,
+                           grid_buf1,grid_buf2,MPI_DOUBLE);
+    else
+      grid3d->reverse_comm(Grid3d::FIX,this,nvalues+1,sizeof(double),0,
+                           grid_buf1,grid_buf2,MPI_DOUBLE);
+  }
+
+  // deallocate local comm buffers b/c new ones will be allocated
 
   memory->destroy(grid_buf1);
   memory->destroy(grid_buf2);
-  memory->destroy3d_offset(T_electron_old, nzlo_out, nylo_out, nxlo_out);
-  memory->destroy3d_offset(net_energy_transfer, nzlo_out, nylo_out, nxlo_out);
 
-  // make copy of ptrs to grid data which does need to persist
+  // make copy of ptrs to grid data which needs to persist
 
-  grid_previous = grid;
+  if (dimension == 2) grid2d_previous = grid2d;
-  T_electron_previous = T_electron;
+  else grid3d_previous = grid3d;
+  
   nxlo_out_previous = nxlo_out;
   nylo_out_previous = nylo_out;
   nzlo_out_previous = nzlo_out;
 
-  // allocate new per-grid data for new decomposition
+  grid_sample_previous = clone_one_grid(grid_sample);
+  grid_nfreq_previous = clone_one_grid(grid_nfreq);
+  if (aveflag == RUNNING || aveflag == WINDOW)
+    grid_running_previous = clone_one_grid(grid_running);
+  if (aveflag == WINDOW) {
+    grid_window_previous = new GridData*[nwindow];
+    for (int i = 0; i < nwindow; i++)
+      grid_window_previous[i] = clone_one_grid(grid_window[i]);
+  }
+
+  // allocate grid instance and grid data for new decomposition
 
   allocate_grid();
 
+  grid_sample = allocate_one_grid();
+  grid_nfreq = allocate_one_grid();
+  if (aveflag == RUNNING || aveflag == WINDOW) grid_nfreq = allocate_one_grid();
+  if (aveflag == WINDOW) {
+    grid_window = new GridData*[nwindow];
+    for (int i = 0; i < nwindow; i++)
+      grid_window[i] = allocate_one_grid();
+  }
+
   // perform remap from previous decomp to new decomp
+  // nper = # of remapped values per grid cell
+  //   depends on atom vs grid mode, and running/window options
 
   int nremap_buf1,nremap_buf2;
-  grid->setup_remap(grid_previous,nremap_buf1,nremap_buf2);
+  if (dimension == 2)
+    grid2d->setup_remap(grid2d_previous,nremap_buf1,nremap_buf2);
+  else
+    grid3d->setup_remap(grid3d_previous,nremap_buf1,nremap_buf2);
 
-  double *remap_buf1,*remap_buf2;
+  int n = 2;                                          // grid_sample & grid_nfreq
-  memory->create(remap_buf1, nremap_buf1, "ave/grid:remap_buf1");
+  if (aveflag == RUNNING || aveflag == WINDOW) n++;   // grid_running
-  memory->create(remap_buf2, nremap_buf2, "ave/grid:remap_buf2");
+  if (aveflag == WINDOW) n += nwindow;                // grid_window
+  int nper = n*nvalues;
+  if (modeatom) nper += n;
 
-  grid->remap(Grid3d::FIX,this,1,sizeof(double),remap_buf1,remap_buf2,MPI_DOUBLE);
+  double *remap_buf1 = nullptr;
+  double *remap_buf2 = nullptr;
+  if (nremap_buf1) memory->create(remap_buf1, nper*nremap_buf1, "ave/grid:remap_buf1");
+  if (nremap_buf2) memory->create(remap_buf2, nper*nremap_buf2, "ave/grid:remap_buf2");
+
+  if (dimension == 2)
+    grid2d->remap(Grid2d::FIX,this,nper,sizeof(double),remap_buf1,remap_buf2,MPI_DOUBLE);
+  else
+    grid3d->remap(Grid3d::FIX,this,nper,sizeof(double),remap_buf1,remap_buf2,MPI_DOUBLE);
 
   memory->destroy(remap_buf1);
   memory->destroy(remap_buf2);
 
-  // delete grid data and grid for previous decomposition
+  // delete grid instance and grid data for previous decomposition
 
-  memory->destroy3d_offset(T_electron_previous,
+  if (dimension == 2) delete grid2d_previous;
-                           nzlo_out_previous, nylo_out_previous,
+  else delete grid3d_previous;
-                           nxlo_out_previous);
-  delete grid_previous;
 
-  // zero new net_energy_transfer
+  deallocate_one_grid(grid_sample_previous,nxlo_out_previous,nylo_out_previous,nzlo_out_previous);
-  // in case compute_vector accesses it on timestep 0
+  deallocate_one_grid(grid_nfreq_previous,nxlo_out_previous,nylo_out_previous,nzlo_out_previous);
+  if (aveflag == RUNNING || aveflag == WINDOW)
+    deallocate_one_grid(grid_running_previous,nxlo_out_previous,nylo_out_previous,nzlo_out_previous);
+  if (aveflag == WINDOW) {
+    for (int i = 0; i < nwindow; i++)
+      deallocate_one_grid(grid_window_previous[i],nxlo_out_previous,nylo_out_previous,nzlo_out_previous);
+    delete [] grid_window_previous;
+  }
 
-  outflag = 0;
+  // set output data in case load balance fix comes after fix ave/grid
-  memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0,
+
-         ngridout*sizeof(double));
+  output_grid(grid_nfreq);
-  */
 }
 
 /* ----------------------------------------------------------------------

src/fix_ave_grid.h

@@ -89,6 +89,14 @@ class FixAveGrid : public Fix {
   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;
+  GridData *grid_sample_previous,*grid_nfreq_previous,*grid_running_previous;
+  GridData **grid_window_previous;
+
   int **bin;
   int *skip;
   int maxatom;
@@ -103,8 +111,14 @@ class FixAveGrid : public Fix {
   void normalize_grid(int, GridData *);
   void normalize_count(int, GridData *);
 
-  void allocate_grid(GridData *);
+  void allocate_grid();
-  void deallocate_grid(GridData *);
+  GridData *allocate_one_grid();
+  GridData *clone_one_grid(GridData *);
+  void deallocate_one_grid(GridData *, int, int, int);
+
+  int pack_one_grid(GridData *, int, double *);
+  int unpack_one_grid(double *, GridData *, int);
+  
   double size_grid(GridData *);
   void zero_grid(GridData *);
   void copy_grid(GridData *, GridData *);

src/grid2d.cpp

@@ -423,14 +423,14 @@ void Grid2d::initialize()
      this includes sub-domain lo boundary but excludes hi boundary
    ngrid = extent of global grid in a dimension
      indices into the global grid range from 0 to Ngrid-1 in that dim
-   shift factor determines position of grid pt within grid cell
+   shift determines position of grid pt within grid cell
      shift = 0.5 for cell center, 0.0 for lower-left corner
    extra = 0 if grid exactly covers the simulation box
    extra = 1 if grid extends beyond the +y boundary by yfactor
      effectively maps proc partitions to the box-size subset of the grid
    lo/hi = inclusive lo/hi bounds for brick of global grid cells I own
-     lo grid index = first grid pt >= fraclo bound
+     lo grid index = first grid pt >= fraclo*bound
-     hi grid index = last grid pt < frachi bound
+     hi grid index = last grid pt < frachi*bound
      if proc owns no grid cells in a dim, then inlo > inhi
    special case: 2 procs share boundary which a grid point is exactly on
      2 if test equalties insure a consistent decision as to which proc owns it
@@ -441,14 +441,14 @@ void Grid2d::partition_grid(int ngrid, double fraclo, double frachi,
 {
   if (extra == 0) {
     lo = static_cast<int> (fraclo * ngrid);
-    while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
+    while (lo+shift < fraclo*ngrid) lo++;
     hi = static_cast<int> (frachi * ngrid);
-    while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
+    while (hi+shift >= frachi*ngrid) hi--;
   } else {
     lo = static_cast<int> (fraclo * ngrid/yfactor);
-    while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
+    while (lo+shift < fraclo*ngrid/yfactor) lo++;
     hi = static_cast<int> (frachi * ngrid/yfactor);
-    while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
+    while (hi+shift >= frachi*ngrid/yfactor) hi--;
   }
 }
 
@@ -556,8 +556,8 @@ void Grid2d::extract_comm_info()
   layout = comm->layout;
 
   // for non TILED layout:
-  // proc xyz lohi = my 64neighbor procs in this MPI_Comm
+  // proc xyz lohi = my 4 neighbor procs in this MPI_Comm
-  // NOTE: will need special logic for MSM case with different MPI_Comm
+  //   these proc IDs can be overridden by caller using set_proc_neighs()
   // xyz split = copy of 1d vectors in Comm
   // grid2proc = copy of 3d array in Comm
 
@@ -609,7 +609,7 @@ void Grid2d::extract_comm_info()
      nbuf2 = larget of sum of all packs or unpacks in Send or Recv
    for brick comm, nbuf1 = nbuf2
    for tiling comm, nbuf2 >= nbuf2
-   nbuf1,nbuf2 are counts of grid points
+   nbuf1,nbuf2 are counts of grid cells
      caller converts them to message sizes for grid data it stores
 ------------------------------------------------------------------------- */
 
@@ -1330,11 +1330,14 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
   int noverlap_old = compute_overlap(0,oldbox,pbc,overlap_old);
 
   // use overlap_old to construct send and copy lists
+  // skip overlaps that contain no grid cells
   
   self_remap = 0;
 
   nsend_remap = 0;
   for (m = 0; m < noverlap_old; m++) {
+    box = overlap_old[m].box;
+    if (box[0] > box[1] || box[2] > box[3]) continue;
     if (overlap_old[m].proc == me) self_remap =1;
     else nsend_remap++;
   }
@@ -1344,6 +1347,7 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
   nsend_remap = 0;
   for (m = 0; m < noverlap_old; m++) {
     box = overlap_old[m].box;
+    if (box[0] > box[1] || box[2] > box[3]) continue;
     if (overlap_old[m].proc == me) {
       copy_remap.npack =
 	old->indices(copy_remap.packlist,box[0],box[1],box[2],box[3]);
@@ -1368,17 +1372,22 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
   int noverlap_new = old->compute_overlap(0,newbox,pbc,overlap_new);
 
   // use overlap_new to construct recv and copy lists
+  // skip overlaps that contain no grid cells
   // set offsets for Recv data
 
   nrecv_remap = 0;
-  for (m = 0; m < noverlap_new; m++)
+  for (m = 0; m < noverlap_new; m++) {
+    box = overlap_new[m].box;
+    if (box[0] > box[1] || box[2] > box[3]) continue;
     if (overlap_new[m].proc != me) nrecv_remap++;
-
+  }
+  
   recv_remap = new Recv[nrecv_remap];
 
   nrecv_remap = 0;
   for (m = 0; m < noverlap_new; m++) {
     box = overlap_new[m].box;
+    if (box[0] > box[1] || box[2] > box[3]) continue;
     if (overlap_new[m].proc == me) {
       copy_remap.nunpack =
 	indices(copy_remap.unpacklist,box[0],box[1],box[2],box[3]);
@@ -1636,14 +1645,23 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
   noverlap_list = maxoverlap_list = 0;
   overlap_list = nullptr;
 
+  // skip overlap check if box contains no grid cells
+
+  if (box[0] > box[1] || box[2] > box[3]) {
+    overlap = overlap_list;
+    return noverlap_list;
+  }
+
+  // test obox against appropriate layout
+  
   if (layout != Comm::LAYOUT_TILED) {
 
     // find comm->procgrid indices in each dim for box bounds
     
-    int iproclo = proc_index_uniform(box[0],nx,0,xsplit);
+    int iproclo = proc_index_uniform(box[0],nx,shift_grid,0,xsplit);
-    int iprochi = proc_index_uniform(box[1],nx,0,xsplit);
+    int iprochi = proc_index_uniform(box[1],nx,shift_grid,0,xsplit);
-    int jproclo = proc_index_uniform(box[2],ny,1,ysplit);
+    int jproclo = proc_index_uniform(box[2],ny,shift_grid,1,ysplit);
-    int jprochi = proc_index_uniform(box[3],ny,1,ysplit);
+    int jprochi = proc_index_uniform(box[3],ny,shift_grid,1,ysplit);
 
     // compute extent of overlap of box with with each proc's obox
 
@@ -1674,6 +1692,9 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
         obox[3] = ny-1;
         
         partition_tiled(overlap_list[m].proc,0,nprocs-1,obox);
+
+        if (me == 1) printf("OBOX: proc %d obox %d %d: %d %d\n",
+                            overlap_list[m].proc,obox[0],obox[1],obox[2],obox[3]);
         
         overlap_list[m].box[0] = MAX(box[0],obox[0]);
         overlap_list[m].box[1] = MIN(box[1],obox[1]);
@@ -1877,29 +1898,33 @@ int Grid2d::indices(int *&list, int xlo, int xhi, int ylo, int yhi)
    find the comm->procgrid index = which proc owns the igrid index
    igrid = grid index (0 to N-1) in dim
    n = # of grid points in dim
+   shift determines position of grid pt within grid cell
+     shift = 0.5 for cell center, 0.0 for lower-left corner
    dim = which dimension (0,1)
    split = comm->x/y/z split for fractional bounds of each proc domain
 ------------------------------------------------------------------------- */
 
-int Grid2d::proc_index_uniform(int igrid, int n, int dim, double *split)
+int Grid2d::proc_index_uniform(int igrid, int n, double shift, int dim, double *split)
 {
-  int gridlo,gridhi;
+  int lo,hi;
   double fraclo,frachi;
 
   // loop over # of procs in this dime
   // compute the grid bounds for that proc
   // if igrid falls within those bounds, return m = proc index
-
+  // same logic as in partition_grid()
+  
   int m;
   for (m = 0; m < comm->procgrid[dim]; m++) {
     fraclo = split[m];
     frachi = split[m+1];
-    gridlo = static_cast<int> (fraclo * n);
+    
-    if (1.0*gridlo != fraclo*n) gridlo++;
+    lo = static_cast<int> (fraclo * n);
-    gridhi = static_cast<int> (frachi * n);
+    while (lo+shift < fraclo*n) lo++;
-    if (1.0*gridhi == frachi*n) gridhi--;
+    hi = static_cast<int> (frachi * n);
+    if (hi+shift >= frachi*n) hi--;
 
-    if (igrid >= gridlo && igrid <= gridhi) break;
+    if (igrid >= lo && igrid <= hi) break;
   }
 
   return m;

src/grid2d.h

@@ -257,7 +257,7 @@ protected:
   void deallocate_remap();
 
   int indices(int *&, int, int, int, int);
-  int proc_index_uniform(int, int, int, double *);
+  int proc_index_uniform(int, int, double, int, double *);
   void partition_tiled(int, int, int, int *);
 };
 

src/grid3d.cpp

@@ -361,7 +361,8 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
                         int &ozlo, int &ozhi)
 {
   // owned grid cells = those whose grid point is within proc subdomain
-  
+  // shift_grid = 0.5 for grid point at cell center, 0.0 for lower-left corner
+
   double fraclo,frachi;
   
   if (comm->layout != Comm::LAYOUT_TILED) {
@@ -469,14 +470,14 @@ void Grid3d::initialize()
      this includes sub-domain lo boundary but excludes hi boundary
    ngrid = extent of global grid in a dimension
      indices into the global grid range from 0 to Ngrid-1 in that dim
-   shift factor determines position of grid pt within grid cell
+   shift determines position of grid pt within grid cell
      shift = 0.5 for cell center, 0.0 for lower-left corner
    extra = 0 if grid exactly covers the simulation box
    extra = 1 if grid extends beyond the +z boundary by zfactor (PPPM slab)
      effectively maps proc partitions to the box-size subset of the grid
    lo/hi = inclusive lo/hi bounds for brick of global grid cells I own
-     lo grid index = first grid pt >= fraclo bound
+     lo grid index = first grid pt >= fraclo*bound
-     hi grid index = last grid pt < frachi bound
+     hi grid index = last grid pt < frachi*bound
      if proc owns no grid cells in a dim, then inlo > inhi
    special case: 2 procs share boundary which a grid point is exactly on
      2 if test equalties insure a consistent decision as to which proc owns it
@@ -487,14 +488,14 @@ void Grid3d::partition_grid(int ngrid, double fraclo, double frachi,
 {
   if (extra == 0) {
     lo = static_cast<int> (fraclo * ngrid);
-    while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
+    while (lo+shift < fraclo*ngrid) lo++;
     hi = static_cast<int> (frachi * ngrid);
-    while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
+    while (hi+shift >= frachi*ngrid) hi--;
   } else {
     lo = static_cast<int> (fraclo * ngrid/zfactor);
-    while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
+    while (lo+shift < fraclo*ngrid/zfactor) lo++;
     hi = static_cast<int> (frachi * ngrid/zfactor);
-    while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
+    while (hi+shift >= frachi*ngrid/zfactor) hi--;
   }
 }
 
@@ -665,7 +666,7 @@ void Grid3d::extract_comm_info()
 
 /* ----------------------------------------------------------------------
    setup commmunication of owned/ghost grid cells
-     either for brick decomp or tiling decomp
+     either for brick decomp or tiled decomp
    return sizes of two buffers needed for communication
      nbuf1 = largest pack or unpack in any Send or Recv or Copy
      nbuf2 = larget of sum of all packs or unpacks in Send or Recv
@@ -1498,11 +1499,14 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
   int noverlap_old = compute_overlap(0,oldbox,pbc,overlap_old);
 
   // use overlap_old to construct send and copy lists
-  
+  // skip overlaps that contain no grid cells
+
   self_remap = 0;
 
   nsend_remap = 0;
   for (m = 0; m < noverlap_old; m++) {
+    box = overlap_old[m].box;
+    if (box[0] > box[1] || box[2] > box[3] || box[4] > box[5]) continue;
     if (overlap_old[m].proc == me) self_remap = 1;
     else nsend_remap++;
   }
@@ -1512,6 +1516,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
   nsend_remap = 0;
   for (m = 0; m < noverlap_old; m++) {
     box = overlap_old[m].box;
+    if (box[0] > box[1] || box[2] > box[3] || box[4] > box[5]) continue;
     if (overlap_old[m].proc == me) {
       copy_remap.npack =
 	old->indices(copy_remap.packlist,
@@ -1537,17 +1542,22 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
   int noverlap_new = old->compute_overlap(0,newbox,pbc,overlap_new);
 
   // use overlap_new to construct recv and copy lists
+  // skip overlaps that contain no grid cells
   // set offsets for Recv data
 
   nrecv_remap = 0;
-  for (m = 0; m < noverlap_new; m++)
+  for (m = 0; m < noverlap_new; m++) {
+    box = overlap_new[m].box;
+    if (box[0] > box[1] || box[2] > box[3] || box[4] > box[5]) continue;
     if (overlap_new[m].proc != me) nrecv_remap++;
-
+  }
+  
   recv_remap = new Recv[nrecv_remap];
 
   nrecv_remap = 0;
   for (m = 0; m < noverlap_new; m++) {
     box = overlap_new[m].box;
+    if (box[0] > box[1] || box[2] > box[3] || box[4] > box[5]) continue;
     if (overlap_new[m].proc == me) {
       copy_remap.nunpack =
 	indices(copy_remap.unpacklist,
@@ -1810,16 +1820,23 @@ int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
   noverlap_list = maxoverlap_list = 0;
   overlap_list = nullptr;
 
+    // skip overlap check if box contains no grid cells
+
+  if (box[0] > box[1] || box[2] > box[3] || box[4] > box[5]) {
+    overlap = overlap_list;
+    return noverlap_list;
+  }
+
   if (layout != Comm::LAYOUT_TILED) {
 
     // find comm->procgrid indices in each dim for box bounds
     
-    int iproclo = proc_index_uniform(box[0],nx,0,xsplit);
+    int iproclo = proc_index_uniform(box[0],nx,shift_grid,0,xsplit);
-    int iprochi = proc_index_uniform(box[1],nx,0,xsplit);
+    int iprochi = proc_index_uniform(box[1],nx,shift_grid,0,xsplit);
-    int jproclo = proc_index_uniform(box[2],ny,1,ysplit);
+    int jproclo = proc_index_uniform(box[2],ny,shift_grid,1,ysplit);
-    int jprochi = proc_index_uniform(box[3],ny,1,ysplit);
+    int jprochi = proc_index_uniform(box[3],ny,shift_grid,1,ysplit);
-    int kproclo = proc_index_uniform(box[4],nz,2,zsplit);
+    int kproclo = proc_index_uniform(box[4],nz,shift_grid,2,zsplit);
-    int kprochi = proc_index_uniform(box[5],nz,2,zsplit);
+    int kprochi = proc_index_uniform(box[5],nz,shift_grid,2,zsplit);
 
     // compute extent of overlap of box with with each proc's obox
     
@@ -2073,30 +2090,34 @@ int Grid3d::indices(int *&list,
 /* ----------------------------------------------------------------------
    find the comm->procgrid index = which proc owns the igrid index
    igrid = grid index (0 to N-1) in dim
-   n = # of grid points in dim
+   n = # of grid points in dim 
+   shift determines position of grid pt within grid cell
+     shift = 0.5 for cell center, 0.0 for lower-left corner
    dim = which dimension (0,1,2)
    split = comm->x/y/z split for fractional bounds of each proc domain
 ------------------------------------------------------------------------- */
 
-int Grid3d::proc_index_uniform(int igrid, int n, int dim, double *split)
+int Grid3d::proc_index_uniform(int igrid, int n, double shift, int dim, double *split)
 {
-  int gridlo,gridhi;
+  int lo,hi;
   double fraclo,frachi;
 
   // loop over # of procs in this dime
   // compute the grid bounds for that proc
   // if igrid falls within those bounds, return m = proc index
+  // same logic as in partition_grid()
 
   int m;
   for (m = 0; m < comm->procgrid[dim]; m++) {
     fraclo = split[m];
     frachi = split[m+1];
-    gridlo = static_cast<int> (fraclo * n);
-    if (1.0*gridlo != fraclo*n) gridlo++;
-    gridhi = static_cast<int> (frachi * n);
-    if (1.0*gridhi == frachi*n) gridhi--;
 
-    if (igrid >= gridlo && igrid <= gridhi) break;
+    lo = static_cast<int> (fraclo * n);
+    while (lo+shift < fraclo*n) lo++;
+    hi = static_cast<int> (frachi * n);
+    if (hi+shift >= frachi*n) hi--;
+
+    if (igrid >= lo && igrid <= hi) break;
   }
 
   return m;
@@ -2108,8 +2129,7 @@ int Grid3d::proc_index_uniform(int igrid, int n, int dim, double *split)
    return box = lo/hi bounds of proc's box in 3 dims
 ------------------------------------------------------------------------- */
 
-void Grid3d::partition_tiled(int proc, int proclower, int procupper,
+void Grid3d::partition_tiled(int proc, int proclower, int procupper, int *box)
-                             int *box)
 {
   // end recursion when partition is a single proc
 

src/grid3d.h

@@ -262,7 +262,7 @@ class Grid3d : protected Pointers {
   void deallocate_remap();
   
   int indices(int *&, int, int, int, int, int, int);
-  int proc_index_uniform(int, int, int, double *);
+  int proc_index_uniform(int, int, double, int, double *);
   void partition_tiled(int, int, int, int *);
 };