From b6583eb681a0930601ad4c9f420dbb0a0a37658a Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 23 Aug 2022 09:58:58 -0600 Subject: [PATCH] debug of fix ave/grid --- src/fix_ave_grid.cpp | 188 ++++++++++++++++++++++++++++++++++++++----- src/fix_ave_grid.h | 6 +- src/grid2d.cpp | 11 ++- src/grid3d.cpp | 11 ++- 4 files changed, 193 insertions(+), 23 deletions(-) diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index 90945433b1..163977c075 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -407,7 +407,12 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : nxlo_in, nxhi_in, nylo_in, nyhi_in, nxlo_out, nxhi_out, nylo_out, nyhi_out); + // ngrid_buf12 converted to nvalues + count + grid2d->setup(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"); @@ -419,7 +424,12 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out); + // ngrid_buf12 converted to nvalues + count + grid3d->setup(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"); @@ -791,7 +801,7 @@ void FixAveGrid::atom2grid() } skip[i] = 0; - count2d[iy][ix] += 1.0; + count2d_sample[iy][ix] += 1.0; bin[i][0] = iy; bin[i][1] = ix; } @@ -824,7 +834,7 @@ void FixAveGrid::atom2grid() } skip[i] = 0; - count3d[iz][iy][ix] += 1.0; + count3d_sample[iz][iy][ix] += 1.0; bin[i][0] = iz; bin[i][1] = iy; bin[i][2] = ix; @@ -857,7 +867,7 @@ void FixAveGrid::atom2grid() } } else for (i = 0; i < nlocal; i++) { - if (skip[i]) + if (!skip[i]) array2d_sample[bin[i][0]][bin[i][1]][m] += attribute[i][j]; } } else { @@ -897,7 +907,7 @@ void FixAveGrid::atom2grid() } } else for (i = 0; i < nlocal; i++) { - if (skip[i]) { + if (!skip[i]) { if (which[m] == ArgInfo::DENSITY_NUMBER) one = 1.0; else if (rmass) one = rmass[i]; else one = mass[type[i]]; @@ -953,7 +963,7 @@ void FixAveGrid::atom2grid() } } else for (i = 0; i < nlocal; i++) { - if (skip[i]) { + if (!skip[i]) { vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; if (rmass) one = rmass[i]; else one = mass[type[i]]; @@ -1292,35 +1302,35 @@ void FixAveGrid::allocate_grid() void FixAveGrid::deallocate_grid() { memory->destroy2d_offset(vec2d,nylo_out,nxlo_out); - memory->destroy2d_offset(array2d,nylo_out,nxlo_out); + memory->destroy3d_offset_last(array2d,nylo_out,nxlo_out); memory->destroy2d_offset(count2d,nylo_out,nxlo_out); memory->destroy2d_offset(vec2d_sample,nylo_out,nxlo_out); - memory->destroy2d_offset(array2d_sample,nylo_out,nxlo_out); + memory->destroy3d_offset_last(array2d_sample,nylo_out,nxlo_out); memory->destroy2d_offset(count2d_sample,nylo_out,nxlo_out); memory->destroy2d_offset(vec2d_epoch,nylo_out,nxlo_out); - memory->destroy2d_offset(array2d_epoch,nylo_out,nxlo_out); + memory->destroy3d_offset_last(array2d_epoch,nylo_out,nxlo_out); memory->destroy2d_offset(count2d_epoch,nylo_out,nxlo_out); memory->destroy2d_offset(vec2d_running,nylo_out,nxlo_out); - memory->destroy2d_offset(array2d_running,nylo_out,nxlo_out); + memory->destroy3d_offset_last(array2d_running,nylo_out,nxlo_out); memory->destroy2d_offset(count2d_running,nylo_out,nxlo_out); memory->destroy3d_offset(vec3d,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(array3d,nzlo_out,nylo_out,nxlo_out); + memory->destroy4d_offset_last(array3d,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(count3d,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(vec3d_sample,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(array3d_sample,nzlo_out,nylo_out,nxlo_out); + memory->destroy4d_offset_last(array3d_sample,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(count3d_sample,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(vec3d_epoch,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(array3d_epoch,nzlo_out,nylo_out,nxlo_out); + memory->destroy4d_offset_last(array3d_epoch,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(count3d_epoch,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(vec3d_running,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(array3d_running,nzlo_out,nylo_out,nxlo_out); + memory->destroy4d_offset_last(array3d_running,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(count3d_running,nzlo_out,nylo_out,nxlo_out); } @@ -1448,6 +1458,150 @@ void FixAveGrid::copy_epoch_to_sample() } } +/* ---------------------------------------------------------------------- + copy sample grid values to output grid, just for owned grid cells + if ATOM mode, also copy per-cell counts +------------------------------------------------------------------------- */ + +void FixAveGrid::copy_sample_to_output() +{ + int ix,iy,iz,m; + + if (!ngridout) return; + + if (dimension == 2) { + if (nvalues == 1) { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec2d[iy][ix] = vec2d_sample[iy][ix]; + } else { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array2d[iy][ix][m] = array2d_sample[iy][ix][m]; + } + if (modeatom) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count2d[iy][ix] = count2d_sample[iy][ix]; + + } else { + if (nvalues == 1) { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec3d[iz][iy][ix] = vec3d_sample[iz][iy][ix]; + } else { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array3d[iz][iy][ix][m] = array3d_sample[iz][iy][ix][m]; + } + if (modeatom) + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count3d[iz][iy][ix] = count3d_sample[iz][iy][ix]; + } +} + +/* ---------------------------------------------------------------------- + copy running grid values to output grid, just for owned grid cells + if ATOM mode, also copy per-cell counts +------------------------------------------------------------------------- */ + +void FixAveGrid::copy_running_to_output() +{ + int ix,iy,iz,m; + + if (!ngridout) return; + + if (dimension == 2) { + if (nvalues == 1) { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec2d[iy][ix] = vec2d_running[iy][ix]; + } else { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array2d[iy][ix][m] = array2d_running[iy][ix][m]; + } + if (modeatom) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count2d[iy][ix] = count2d_running[iy][ix]; + + } else { + if (nvalues == 1) { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec3d[iz][iy][ix] = vec3d_running[iz][iy][ix]; + } else { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array3d[iz][iy][ix][m] = array3d_running[iz][iy][ix][m]; + } + if (modeatom) + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count3d[iz][iy][ix] = count3d_running[iz][iy][ix]; + } +} + +/* ---------------------------------------------------------------------- + sum sample grid values to running grid, just for owned grid cells + if ATOM mode, also copy per-cell counts +------------------------------------------------------------------------- */ + +void FixAveGrid::sum_sample_to_running() +{ + int ix,iy,iz,m; + + if (!ngridout) return; + + if (dimension == 2) { + if (nvalues == 1) { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec2d_running[iy][ix] += vec2d_sample[iy][ix]; + } else { + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array2d_running[iy][ix][m] += array2d_sample[iy][ix][m]; + } + if (modeatom) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count2d_running[iy][ix] += count2d_sample[iy][ix]; + + } else { + if (nvalues == 1) { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + vec3d_running[iz][iy][ix] += vec3d_sample[iz][iy][ix]; + } else { + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + for (m = 0; m <= nvalues; m++) + array3d_running[iz][iy][ix][m] += array3d_sample[iz][iy][ix][m]; + } + if (modeatom) + for (iz = nylo_in; iz <= nyhi_in; iz++) + for (iy = nylo_in; iy <= nyhi_in; iy++) + for (ix = nxlo_in; ix <= nxhi_in; ix++) + count3d_running[iz][iy][ix] += count3d_sample[iz][iy][ix]; + } +} + /* ---------------------------------------------------------------------- normalize grid values and per-cell counts for ATOM mode for owned cells numsamples = normalization factor @@ -1496,7 +1650,7 @@ void FixAveGrid::normalize_atom(int numsamples) norm = mvv2e /((repeat*cdof + adof*count) * boltz); else if (normflag == NONORM) norm = 1.0/invrepeat; - else if (normflag == NONORM) + else norm = 1.0/count; vec2d_sample[iy][ix] *= norm; count2d_sample[iy][ix] *= invrepeat; @@ -1507,7 +1661,7 @@ void FixAveGrid::normalize_atom(int numsamples) for (ix = nxlo_in; ix <= nxhi_in; ix++) { count = count2d_sample[iy][ix]; if (count) { - for (m = 0; m <= nvalues; m++) { + for (m = 0; m < nvalues; m++) { if (which[m] == ArgInfo::DENSITY_NUMBER) norm = 1.0 / (binvol * repeat); else if (which[m] == ArgInfo::DENSITY_MASS) @@ -1551,7 +1705,7 @@ void FixAveGrid::normalize_atom(int numsamples) for (ix = nxlo_in; ix <= nxhi_in; ix++) { count = count3d_sample[iz][iy][ix]; if (count) { - for (m = 0; m <= nvalues; m++) { + for (m = 0; m < nvalues; m++) { if (which[m] == ArgInfo::DENSITY_NUMBER) norm = 1.0 / (binvol * repeat); else if (which[m] == ArgInfo::DENSITY_MASS) @@ -1641,8 +1795,6 @@ void FixAveGrid::pack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *lis } } else { for (i = 0; i < nlist; i++) { - printf("PACK i %d nlist %d, listI %d nvalues %d count %d\n", - i,nlist,list[i],nvalues,count[list[i]]); buf[m++] = count[list[i]]; values = &data[nvalues*list[i]]; for (j = 0; j < nvalues; j++) diff --git a/src/fix_ave_grid.h b/src/fix_ave_grid.h index e8efee4eed..cce7c5bb52 100644 --- a/src/fix_ave_grid.h +++ b/src/fix_ave_grid.h @@ -113,9 +113,9 @@ class FixAveGrid : public Fix { double ***, double ***, double ****); void sum_sample_to_epoch(); void copy_epoch_to_sample(); - void sum_sample_to_running() {} - void copy_sample_to_output() {} - void copy_running_to_output() {} + void sum_sample_to_running(); + void copy_sample_to_output(); + void copy_running_to_output(); void copy_sample_to_window(int) {} void subtract_window_from_running() {} diff --git a/src/grid2d.cpp b/src/grid2d.cpp index 8aef2acbe4..9077493269 100644 --- a/src/grid2d.cpp +++ b/src/grid2d.cpp @@ -366,7 +366,16 @@ void Grid2d::get_box(int dim, double &lo, double &delta) delta = prd[dim] / ngrid[dim]; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + return sizes of two buffers needed for communication + either on regular grid or procs or irregular tiling + 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 + for regular comm, nbuf1 = nbuf2 + for irregular comm, nbuf2 >= nbuf2 + nbuf1,nbuf2 are just count of grid points + caller converts them to message size for grid data it stores +------------------------------------------------------------------------- */ void Grid2d::setup(int &nbuf1, int &nbuf2) { diff --git a/src/grid3d.cpp b/src/grid3d.cpp index fb7334e1f1..51d2ebd247 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -399,7 +399,16 @@ void Grid3d::get_box(int dim, double &lo, double &delta) delta = prd[dim] / ngrid[dim]; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + return sizes of two buffers needed for communication + either on regular grid or procs or irregular tiling + 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 + for regular comm, nbuf1 = nbuf2 + for irregular comm, nbuf2 >= nbuf2 + nbuf1,nbuf2 are just count of grid points + caller converts them to message size for grid data it stores +------------------------------------------------------------------------- */ void Grid3d::setup(int &nbuf1, int &nbuf2) {