From 3dab9cf8d86a72953385672c9b95aaca0c5efa7d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 4 Nov 2022 21:41:58 -0600 Subject: [PATCH] bug fixes to grid remapping --- src/dump_grid.cpp | 3 + src/fix_ave_grid.cpp | 444 ++++++++++++++++++++++++++++++------------- src/fix_ave_grid.h | 18 +- src/grid2d.cpp | 73 ++++--- src/grid2d.h | 2 +- src/grid3d.cpp | 76 +++++--- src/grid3d.h | 2 +- 7 files changed, 425 insertions(+), 193 deletions(-) diff --git a/src/dump_grid.cpp b/src/dump_grid.cpp index 3a32b95bd4..aa4d37c1e6 100644 --- a/src/dump_grid.cpp +++ b/src/dump_grid.cpp @@ -26,6 +26,9 @@ #include "region.h" #include "update.h" +// DEBUG +#include "comm.h" + #include using namespace LAMMPS_NS; diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index dcc1016729..16f62d57bc 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -29,9 +29,6 @@ #include "update.h" #include "variable.h" -// DEBUG -#include "comm.h" - #include 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; - else if (modegrid) maxdist = 0.0; + allocate_grid(); - 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; - - 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(); + 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] = new GridData(); - } else grid_window = nullptr; + grid_window[i] = allocate_one_grid(); + } - allocate_grid(grid_sample); - allocate_grid(grid_nfreq); - if (aveflag == RUNNING || aveflag == WINDOW) allocate_grid(grid_running); - if (aveflag == WINDOW) - for (int i = 0; i < nwindow; i++) - allocate_grid(grid_window[i]); + // output may occur via dump on timestep 0 + + grid_output = new GridData(); + output_grid(grid_nfreq); // 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; - delete grid3d; + // deallocate Grid class and buffers + + 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_grid(grid_nfreq); - if (aveflag == RUNNING || aveflag == WINDOW) deallocate_grid(grid_running); - if (aveflag == WINDOW) - for (int i = 0; i < nwindow; i++) { - deallocate_grid(grid_window[i]); - delete grid_window[i]; - } + deallocate_one_grid(grid_sample,nxlo_out,nylo_out,nzlo_out); + deallocate_one_grid(grid_nfreq,nxlo_out,nylo_out,nzlo_out); + if (aveflag == RUNNING || aveflag == WINDOW) + deallocate_one_grid(grid_running,nxlo_out,nylo_out,nzlo_out); + if (aveflag == WINDOW) { + for (int i = 0; i < nwindow; 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; - double *src; - //double *src = - // &T_electron_previous[nzlo_out_previous][nylo_out_previous][nxlo_out_previous]; + int i,j,m,iwindow; - 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; - double *dest; - //double *dest = &T_electron[nzlo_out][nylo_out][nxlo_out]; + int i,j,m,iwindow; - 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 - if (comm->me == 0) printf("Remapping grid on step %ld\n",update->ntimestep); + // for ATOM mode, perform ghost to owned grid comm for grid_sample + // 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 - /* - // delete grid data which doesn't need to persist from previous to new decomp + if (modeatom) { + 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; - T_electron_previous = T_electron; + if (dimension == 2) grid2d_previous = grid2d; + 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; - memory->create(remap_buf1, nremap_buf1, "ave/grid:remap_buf1"); - memory->create(remap_buf2, nremap_buf2, "ave/grid:remap_buf2"); + int n = 2; // grid_sample & grid_nfreq + if (aveflag == RUNNING || aveflag == WINDOW) n++; // grid_running + 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, - nzlo_out_previous, nylo_out_previous, - nxlo_out_previous); - delete grid_previous; + if (dimension == 2) delete grid2d_previous; + else delete grid3d_previous; - // zero new net_energy_transfer - // in case compute_vector accesses it on timestep 0 + deallocate_one_grid(grid_sample_previous,nxlo_out_previous,nylo_out_previous,nzlo_out_previous); + 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; - memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0, - ngridout*sizeof(double)); - */ + // set output data in case load balance fix comes after fix ave/grid + + output_grid(grid_nfreq); } /* ---------------------------------------------------------------------- diff --git a/src/fix_ave_grid.h b/src/fix_ave_grid.h index 71f1a34d8a..53c90c2458 100644 --- a/src/fix_ave_grid.h +++ b/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 deallocate_grid(GridData *); + void allocate_grid(); + 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 *); diff --git a/src/grid2d.cpp b/src/grid2d.cpp index 77ae289d63..61ced82c59 100644 --- a/src/grid2d.cpp +++ b/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 - hi grid index = last grid pt < frachi bound + lo grid index = first grid pt >= fraclo*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 (fraclo * ngrid); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + while (lo+shift < fraclo*ngrid) lo++; hi = static_cast (frachi * ngrid); - while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + while (hi+shift >= frachi*ngrid) hi--; } else { lo = static_cast (fraclo * ngrid/yfactor); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + while (lo+shift < fraclo*ngrid/yfactor) lo++; hi = static_cast (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 - // NOTE: will need special logic for MSM case with different MPI_Comm + // proc xyz lohi = my 4 neighbor procs in this 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 iprochi = proc_index_uniform(box[1],nx,0,xsplit); - int jproclo = proc_index_uniform(box[2],ny,1,ysplit); - int jprochi = proc_index_uniform(box[3],ny,1,ysplit); + int iproclo = proc_index_uniform(box[0],nx,shift_grid,0,xsplit); + int iprochi = proc_index_uniform(box[1],nx,shift_grid,0,xsplit); + int jproclo = proc_index_uniform(box[2],ny,shift_grid,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 (fraclo * n); - if (1.0*gridlo != fraclo*n) gridlo++; - gridhi = static_cast (frachi * n); - if (1.0*gridhi == frachi*n) gridhi--; + + lo = static_cast (fraclo * n); + while (lo+shift < fraclo*n) lo++; + hi = static_cast (frachi * n); + if (hi+shift >= frachi*n) hi--; - if (igrid >= gridlo && igrid <= gridhi) break; + if (igrid >= lo && igrid <= hi) break; } return m; diff --git a/src/grid2d.h b/src/grid2d.h index 3ad28b3a4f..3772028f52 100644 --- a/src/grid2d.h +++ b/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 *); }; diff --git a/src/grid3d.cpp b/src/grid3d.cpp index 30056d27c9..ddf2667ad1 100644 --- a/src/grid3d.cpp +++ b/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 - hi grid index = last grid pt < frachi bound + lo grid index = first grid pt >= fraclo*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 (fraclo * ngrid); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + while (lo+shift < fraclo*ngrid) lo++; hi = static_cast (frachi * ngrid); - while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + while (hi+shift >= frachi*ngrid) hi--; } else { lo = static_cast (fraclo * ngrid/zfactor); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + while (lo+shift < fraclo*ngrid/zfactor) lo++; hi = static_cast (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 iprochi = proc_index_uniform(box[1],nx,0,xsplit); - int jproclo = proc_index_uniform(box[2],ny,1,ysplit); - int jprochi = proc_index_uniform(box[3],ny,1,ysplit); - int kproclo = proc_index_uniform(box[4],nz,2,zsplit); - int kprochi = proc_index_uniform(box[5],nz,2,zsplit); + int iproclo = proc_index_uniform(box[0],nx,shift_grid,0,xsplit); + int iprochi = proc_index_uniform(box[1],nx,shift_grid,0,xsplit); + int jproclo = proc_index_uniform(box[2],ny,shift_grid,1,ysplit); + int jprochi = proc_index_uniform(box[3],ny,shift_grid,1,ysplit); + int kproclo = proc_index_uniform(box[4],nz,shift_grid,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 (fraclo * n); - if (1.0*gridlo != fraclo*n) gridlo++; - gridhi = static_cast (frachi * n); - if (1.0*gridhi == frachi*n) gridhi--; - if (igrid >= gridlo && igrid <= gridhi) break; + lo = static_cast (fraclo * n); + while (lo+shift < fraclo*n) lo++; + hi = static_cast (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, - int *box) +void Grid3d::partition_tiled(int proc, int proclower, int procupper, int *box) { // end recursion when partition is a single proc diff --git a/src/grid3d.h b/src/grid3d.h index 73f6bcb5c5..721cdd93f0 100644 --- a/src/grid3d.h +++ b/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 *); };