From 4c29457351b0c2d0198caabfe66a685db3cb717b Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 1 Nov 2022 15:53:39 -0600 Subject: [PATCH] more classes use Grid3d --- src/AMOEBA/amoeba_convolution.cpp | 114 +++------- src/AMOEBA/amoeba_kspace.cpp | 4 + src/EXTRA-FIX/fix_ttm_grid.cpp | 17 +- src/KSPACE/msm.cpp | 4 +- src/KSPACE/pppm.cpp | 166 ++++++-------- src/KSPACE/pppm.h | 2 +- src/KSPACE/pppm_dipole.cpp | 70 ++++-- src/KSPACE/pppm_dipole_spin.cpp | 24 +- src/KSPACE/pppm_disp.cpp | 12 +- src/compute_property_grid.cpp | 6 +- src/dump_grid.cpp | 2 +- src/fix_ave_grid.cpp | 8 +- src/grid3d.cpp | 367 ++++++++++++++++++------------ src/grid3d.h | 37 ++- 14 files changed, 440 insertions(+), 393 deletions(-) diff --git a/src/AMOEBA/amoeba_convolution.cpp b/src/AMOEBA/amoeba_convolution.cpp index eca1c88c4d..ac6a8d161b 100644 --- a/src/AMOEBA/amoeba_convolution.cpp +++ b/src/AMOEBA/amoeba_convolution.cpp @@ -104,82 +104,29 @@ void AmoebaConvolution::reset_grid() void AmoebaConvolution::allocate_grid() { - // global indices of grid range from 0 to N-1 - // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of - // global grid that I own without ghost cells - // both non-tiled and tiled proc layouts use 0-1 fractional subdomain info - - if (comm->layout != Comm::LAYOUT_TILED) { - nxlo_in = static_cast (comm->xsplit[comm->myloc[0]] * nx); - nxhi_in = static_cast (comm->xsplit[comm->myloc[0]+1] * nx) - 1; - nylo_in = static_cast (comm->ysplit[comm->myloc[1]] * ny); - nyhi_in = static_cast (comm->ysplit[comm->myloc[1]+1] * ny) - 1; - nzlo_in = static_cast (comm->zsplit[comm->myloc[2]] * nz); - nzhi_in = static_cast (comm->zsplit[comm->myloc[2]+1] * nz) - 1; - - } else { - nxlo_in = static_cast (comm->mysplit[0][0] * nx); - nxhi_in = static_cast (comm->mysplit[0][1] * nx) - 1; - nylo_in = static_cast (comm->mysplit[1][0] * ny); - nyhi_in = static_cast (comm->mysplit[1][1] * ny) - 1; - nzlo_in = static_cast (comm->mysplit[2][0] * nz); - nzhi_in = static_cast (comm->mysplit[2][1] * nz) - 1; - } - + // maxdist = max distance outside of subbox my owned atom may be // nlower,nupper = stencil size for mapping particles to FFT grid + double maxdist = 0.5*neighbor->skin; int nlower = -(order-1)/2; int nupper = order/2; - // nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of - // global grid that my particles can contribute charge to - // effectively nlo_in,nhi_in + ghost cells - // nlo,nhi = global coords of grid pt to "lower left" of smallest/largest - // position a particle in my box can be at - // dist[3] = particle position bound = subbox + skin/2.0 - // convert to triclinic if necessary - // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping + // Grid3d determines my owned + ghost grid cells + // ghost cell extent depends on maxdist, nlower, nupper + + gc = new Grid3d(lmp,world,nx,ny,nz); + gc->set_distance(maxdist); + gc->set_stencil_atom(-nlower,nupper); - double *prd,*boxlo,*sublo,*subhi; - int triclinic = domain->triclinic; + gc->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); - if (triclinic == 0) { - prd = domain->prd; - boxlo = domain->boxlo; - sublo = domain->sublo; - subhi = domain->subhi; - } else { - prd = domain->prd_lamda; - boxlo = domain->boxlo_lamda; - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } + int nqty = flag3d ? 1 : 2; + int ngc_buf1,ngc_buf2; - double xprd = prd[0]; - double yprd = prd[1]; - double zprd = prd[2]; - - double dist[3] = {0.0,0.0,0.0}; - double cuthalf = 0.5*neighbor->skin; - if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf; - else MathExtra::tribbox(domain->h,cuthalf,&dist[0]); - - int nlo,nhi; - - nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * nx/xprd); - nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * nx/xprd); - nxlo_out = nlo + nlower; - nxhi_out = nhi + nupper; - - nlo = static_cast ((sublo[1]-dist[1]-boxlo[1]) * ny/yprd); - nhi = static_cast ((subhi[1]+dist[1]-boxlo[1]) * ny/yprd); - nylo_out = nlo + nlower; - nyhi_out = nhi + nupper; - - nlo = static_cast ((sublo[2]-dist[2]-boxlo[2]) * nz/zprd); - nhi = static_cast ((subhi[2]+dist[2]-boxlo[2]) * nz/zprd); - nzlo_out = nlo + nlower; - nzhi_out = nhi + nupper; + gc->setup_comm(ngc_buf1,ngc_buf2); + memory->create(gc_buf1,nqty*ngc_buf1,"amoeba:gc_buf1"); + memory->create(gc_buf2,nqty*ngc_buf2,"amoeba:gc_buf2"); // x-pencil decomposition of FFT mesh // global indices range from 0 to N-1 @@ -210,7 +157,7 @@ void AmoebaConvolution::allocate_grid() nzlo_fft = me_z*nz/npez_fft; nzhi_fft = (me_z+1)*nz/npez_fft - 1; - // grid sizes + // FFT grid sizes // nbrick_owned = owned grid points in brick decomp // nbrick_ghosts = owned + ghost grid points in grid decomp // nfft_owned = owned grid points in FFT decomp @@ -226,7 +173,7 @@ void AmoebaConvolution::allocate_grid() ngrid_either = MAX(nbrick_owned,nfft_owned); - // instantiate FFT, Grid3d, and Remap + // instantiate FFT and Remap int tmp; @@ -240,11 +187,6 @@ void AmoebaConvolution::allocate_grid() nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, 0,0,&tmp,0); - gc = new Grid3d(lmp,world,nx,ny,nz, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); - - int nqty = flag3d ? 1 : 2; remap = new Remap(lmp,world, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, @@ -266,12 +208,6 @@ void AmoebaConvolution::allocate_grid() memory->create(grid_fft,ngrid_either,"amoeba:grid_fft"); memory->create(cfft,2*ngrid_either,"amoeba:cfft"); - - int ngc_buf1,ngc_buf2; - gc->setup(ngc_buf1,ngc_buf2); - memory->create(gc_buf1,nqty*ngc_buf1,"amoeba:gc_buf1"); - memory->create(gc_buf2,nqty*ngc_buf2,"amoeba:gc_buf2"); - memory->create(remap_buf,nqty*nfft_owned,"amoeba:remap_buf"); } @@ -279,18 +215,20 @@ void AmoebaConvolution::allocate_grid() void AmoebaConvolution::deallocate_grid() { + delete gc; + + memory->destroy(gc_buf1); + memory->destroy(gc_buf2); + + delete fft1; + delete fft2; + delete remap; + memory->destroy3d_offset(grid_brick,nzlo_out,nylo_out,nxlo_out); memory->destroy4d_offset_last(cgrid_brick,nzlo_out,nylo_out,nxlo_out); memory->destroy(grid_fft); memory->destroy(cfft); - memory->destroy(gc_buf1); - memory->destroy(gc_buf2); memory->destroy(remap_buf); - - delete fft1; - delete fft2; - delete gc; - delete remap; } /* ---------------------------------------------------------------------- diff --git a/src/AMOEBA/amoeba_kspace.cpp b/src/AMOEBA/amoeba_kspace.cpp index 51b206b6d8..c4fb5136b0 100644 --- a/src/AMOEBA/amoeba_kspace.cpp +++ b/src/AMOEBA/amoeba_kspace.cpp @@ -217,6 +217,10 @@ void PairAmoeba::bspline_fill() // NOTE: could subtract off nlpts to start with // NOTE: this is place to check that stencil size does not // go out of bounds relative to igrid for a proc's sub-domain + // similar to PPPM::particle_map() + // subtracting eps is strange, and could mess up the check + // better to check here than in methods like grid_mpole() + // but check needs to be valid for all KSpace terms // NOTE: could convert x -> lamda for entire set of Nlocal atoms domain->x2lamda(x[i],lamda); diff --git a/src/EXTRA-FIX/fix_ttm_grid.cpp b/src/EXTRA-FIX/fix_ttm_grid.cpp index b5b57ca5e1..08f1e43c32 100644 --- a/src/EXTRA-FIX/fix_ttm_grid.cpp +++ b/src/EXTRA-FIX/fix_ttm_grid.cpp @@ -485,10 +485,11 @@ void FixTTMGrid::reset_grid() int tmp[12]; double maxdist = 0.5 * neighbor->skin; - Grid3d *gridnew = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, - maxdist, 1, 0.5, - tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],tmp[5], - tmp[6],tmp[7],tmp[8],tmp[9],tmp[10],tmp[11]); + Grid3d *gridnew = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid); + gridnew->set_distance(maxdist); + gridnew->set_stencil_grid(1,1); + gridnew->setup_grid(tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],tmp[5], + tmp[6],tmp[7],tmp[8],tmp[9],tmp[10],tmp[11]); if (grid->identical(gridnew)) { delete gridnew; @@ -632,9 +633,11 @@ void FixTTMGrid::allocate_grid() { double maxdist = 0.5 * neighbor->skin; - grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 1, 0.5, - nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in, - nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out); + grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid); + grid->set_distance(maxdist); + grid->set_stencil_grid(1,1); + grid->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); ngridown = (nxhi_in - nxlo_in + 1) * (nyhi_in - nylo_in + 1) * (nzhi_in - nzlo_in + 1); diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index a3ef458770..ebc5858f10 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -604,7 +604,7 @@ void MSM::allocate() nxlo_out[0],nxhi_out[0],nylo_out[0], nyhi_out[0],nzlo_out[0],nzhi_out[0]); - gcall->setup(ngcall_buf1,ngcall_buf2); + gcall->setup_comm(ngcall_buf1,ngcall_buf2); npergrid = 1; memory->destroy(gcall_buf1); memory->destroy(gcall_buf2); @@ -636,7 +636,7 @@ void MSM::allocate() procneigh[0][0],procneigh[0][1],procneigh[1][0], procneigh[1][1],procneigh[2][0],procneigh[2][1]); - gc[n]->setup(ngc_buf1[n],ngc_buf2[n]); + gc[n]->setup_comm(ngc_buf1[n],ngc_buf2[n]); npergrid = 1; memory->destroy(gc_buf1[n]); memory->destroy(gc_buf2[n]); diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index f946d44d29..2ad3b15c72 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -292,7 +292,7 @@ void PPPM::init() // or overlap is allowed, then done // else reduce order and try again - Grid3d *gctmp = nullptr; + gc = nullptr; int iteration = 0; while (order >= minorder) { @@ -305,23 +305,28 @@ void PPPM::init() set_grid_local(); if (overlap_allowed) break; - gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); + gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm); + gc->set_distance(0.5*neighbor->skin + qdist); + gc->set_stencil_atom(-nlower,nupper); + gc->set_shift_atom(shiftatom); + gc->set_zfactor(slab_volfactor); + + gc->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); int tmp1,tmp2; - gctmp->setup(tmp1,tmp2); - if (gctmp->ghost_adjacent()) break; - delete gctmp; + gc->setup_comm(tmp1,tmp2); + if (gc->ghost_adjacent()) break; + delete gc; order--; iteration++; } if (order < minorder) error->all(FLERR,"PPPM order < minimum allowed order"); - if (!overlap_allowed && !gctmp->ghost_adjacent()) + if (!overlap_allowed && !gc->ghost_adjacent()) error->all(FLERR,"PPPM grid stencil extends beyond nearest neighbor processor"); - if (gctmp) delete gctmp; + if (gc) delete gc; // adjust g_ewald @@ -739,6 +744,29 @@ void PPPM::compute(int eflag, int vflag) void PPPM::allocate() { + // create ghost grid object for rho and electric field communication + // returns local owned and ghost grid bounds + // setup communication patterns and buffers + + gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm); + gc->set_distance(0.5*neighbor->skin + qdist); + gc->set_stencil_atom(-nlower,nupper); + gc->set_shift_atom(shiftatom); + gc->set_zfactor(slab_volfactor); + + gc->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); + + gc->setup_comm(ngc_buf1,ngc_buf2); + + if (differentiation_flag) npergrid = 1; + else npergrid = 3; + + memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1"); + memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2"); + + // allocate distributed grid data + memory->create3d_offset(density_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:density_brick"); @@ -809,21 +837,6 @@ void PPPM::allocate() nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, 1,0,0,FFT_PRECISION,collective_flag); - - // create ghost grid object for rho and electric field communication - // also create 2 bufs for ghost grid cell comm, passed to Grid3d methods - - gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); - - gc->setup(ngc_buf1,ngc_buf2); - - if (differentiation_flag) npergrid = 1; - else npergrid = 3; - - memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1"); - memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2"); } /* ---------------------------------------------------------------------- @@ -832,6 +845,10 @@ void PPPM::allocate() void PPPM::deallocate() { + delete gc; + memory->destroy(gc_buf1); + memory->destroy(gc_buf2); + memory->destroy3d_offset(density_brick,nzlo_out,nylo_out,nxlo_out); if (differentiation_flag == 1) { @@ -874,9 +891,6 @@ void PPPM::deallocate() delete fft1; delete fft2; delete remap; - delete gc; - memory->destroy(gc_buf1); - memory->destroy(gc_buf2); } /* ---------------------------------------------------------------------- @@ -1309,85 +1323,41 @@ double PPPM::final_accuracy() } /* ---------------------------------------------------------------------- - set local subset of PPPM/FFT grid that I own - n xyz lo/hi in = 3d brick that I own (inclusive) - n xyz lo/hi out = 3d brick + ghost cells in 6 directions (inclusive) - n xyz lo/hi fft = FFT columns that I own (all of x dim, 2d decomp in yz) + set params which determine which owned and ghost cells this proc owns + Grid3d will use these params to partition grid + also partition FFT grid + n xyz lo/hi fft = FFT columns that I own (all of x dim, 2d decomp in yz) ------------------------------------------------------------------------- */ void PPPM::set_grid_local() { - // partition global grid across procs - // n xyz lo/hi in = lower/upper bounds of global grid this proc owns - // indices range from 0 to N-1 inclusive in each dim + // shift values for particle <-> grid mapping + // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 - comm->partition_grid(nx_pppm,ny_pppm,nz_pppm,slab_volfactor, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in); + if (order % 2) shiftatom = 0.5; + else shiftatom = 0.0; + shift = OFFSET + shiftatom; + + if (order % 2) shiftone = 0.0; + else shiftone = 0.5; // nlower,nupper = stencil size for mapping particles to PPPM grid nlower = -(order-1)/2; nupper = order/2; - // shift values for particle <-> grid mapping - // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 - if (order % 2) shift = OFFSET + 0.5; - else shift = OFFSET; - if (order % 2) shiftone = 0.0; - else shiftone = 0.5; - // nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of - // global PPPM grid that my particles can contribute charge to - // effectively nlo_in,nhi_in + ghost cells - // nlo,nhi = index of global grid pt to "lower left" of smallest/largest - // position a particle in my box can be at - // dist[3] = max particle position outside subbox = skin/2.0 + qdist - // qdist = offset due to TIP4P fictitious charge - // convert to triclinic if necessary - // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping - // for slab PPPM, assign z grid as if it were not extended + + // NOTE: still to do: stagger and zperiod effects - double *prd,*sublo,*subhi; + /* - if (triclinic == 0) { - prd = domain->prd; - boxlo = domain->boxlo; - sublo = domain->sublo; - subhi = domain->subhi; - } else { - prd = domain->prd_lamda; - boxlo = domain->boxlo_lamda; - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } - - double xprd = prd[0]; - double yprd = prd[1]; - double zprd = prd[2]; + // extent of zprd when 2d slab mode is selected + double zprd_slab = zprd*slab_volfactor; - double dist[3] = {0.0,0.0,0.0}; - double cuthalf = 0.5*neighbor->skin + qdist; - if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf; - else MathExtra::tribbox(domain->h,cuthalf,&dist[0]); - - int nlo,nhi; - nlo = nhi = 0; - - nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * - nx_pppm/xprd + shift) - OFFSET; - nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * - nx_pppm/xprd + shift) - OFFSET; - nxlo_out = nlo + nlower; - nxhi_out = nhi + nupper; - - nlo = static_cast ((sublo[1]-dist[1]-boxlo[1]) * - ny_pppm/yprd + shift) - OFFSET; - nhi = static_cast ((subhi[1]+dist[1]-boxlo[1]) * - ny_pppm/yprd + shift) - OFFSET; - nylo_out = nlo + nlower; - nyhi_out = nhi + nupper; + // for slab PPPM, assign z grid as if it were not extended nlo = static_cast ((sublo[2]-dist[2]-boxlo[2]) * nz_pppm/zprd_slab + shift) - OFFSET; @@ -1396,6 +1366,7 @@ void PPPM::set_grid_local() nzlo_out = nlo + nlower; nzhi_out = nhi + nupper; + if (stagger_flag) { nxhi_out++; nyhi_out++; @@ -1420,6 +1391,11 @@ void PPPM::set_grid_local() nzhi_out = MIN(nzhi_out,nz_pppm-1); } + */ + + + + // x-pencil decomposition of FFT mesh // global indices range from 0 to N-1 // each proc owns entire x-dimension, clumps of columns in y,z dimensions @@ -1446,18 +1422,22 @@ void PPPM::set_grid_local() nzlo_fft = me_z*nz_pppm/npez_fft; nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1; + // nfft = FFT points in x-pencil FFT decomposition on this proc + + nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) * + (nzhi_fft-nzlo_fft+1); + + + // ngrid = count of PPPM grid pts owned by this proc, including ghosts ngrid = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * (nzhi_out-nzlo_out+1); // count of FFT grids pts owned by this proc, without ghosts - // nfft = FFT points in x-pencil FFT decomposition on this proc // nfft_brick = FFT points in 3d brick-decomposition on this proc // nfft_both = greater of 2 values - nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) * - (nzhi_fft-nzlo_fft+1); int nfft_brick = (nxhi_in-nxlo_in+1) * (nyhi_in-nylo_in+1) * (nzhi_in-nzlo_in+1); nfft_both = MAX(nfft,nfft_brick); @@ -1908,7 +1888,7 @@ void PPPM::make_rho() // loop over my charges, add their contribution to nearby grid points // (nx,ny,nz) = global coords of grid pt to "lower left" of charge // (dx,dy,dz) = distance to "lower left" grid pt - // (mx,my,mz) = global coords of moving stencil pt + // (mx,my,mz) = global indices of moving stencil pt double *q = atom->q; double **x = atom->x; diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index da1278ff4f..6166961cba 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -48,7 +48,7 @@ class PPPM : public KSpace { double volume; double delxinv, delyinv, delzinv, delvolinv; double h_x, h_y, h_z; - double shift, shiftone; + double shift, shiftone, shiftatom; int peratom_allocate_flag; int nxlo_in, nylo_in, nzlo_in, nxhi_in, nyhi_in, nzhi_in; diff --git a/src/KSPACE/pppm_dipole.cpp b/src/KSPACE/pppm_dipole.cpp index d6931d4c0e..a2abc98961 100644 --- a/src/KSPACE/pppm_dipole.cpp +++ b/src/KSPACE/pppm_dipole.cpp @@ -28,6 +28,7 @@ #include "math_const.h" #include "math_special.h" #include "memory.h" +#include "neighbor.h" #include "pair.h" #include "remap_wrap.h" #include "update.h" @@ -188,7 +189,7 @@ void PPPMDipole::init() // or overlap is allowed, then done // else reduce order and try again - Grid3d *gctmp = nullptr; + gc_dipole = nullptr; int iteration = 0; while (order >= minorder) { @@ -201,23 +202,28 @@ void PPPMDipole::init() set_grid_local(); if (overlap_allowed) break; - gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); + gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm); + gc_dipole->set_distance(0.5*neighbor->skin + qdist); + gc_dipole->set_stencil_atom(-nlower,nupper); + gc_dipole->set_shift_atom(shiftatom); + gc_dipole->set_zfactor(slab_volfactor); + + gc_dipole->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); int tmp1,tmp2; - gctmp->setup(tmp1,tmp2); - if (gctmp->ghost_adjacent()) break; - delete gctmp; + gc_dipole->setup_comm(tmp1,tmp2); + if (gc_dipole->ghost_adjacent()) break; + delete gc_dipole; order--; iteration++; } if (order < minorder) error->all(FLERR,"PPPMDipole order < minimum allowed order"); - if (!overlap_allowed && !gctmp->ghost_adjacent()) + if (!overlap_allowed && !gc_dipole->ghost_adjacent()) error->all(FLERR,"PPPMDipole grid stencil extends beyond nearest neighbor processor"); - if (gctmp) delete gctmp; + if (gc_dipole) delete gc_dipole; // adjust g_ewald @@ -436,7 +442,7 @@ void PPPMDipole::compute(int eflag, int vflag) particle_map(); make_rho_dipole(); - + // all procs communicate density values from their ghost cells // to fully sum contribution in their 3d bricks // remap from 3d decomposition to FFT decomposition @@ -528,6 +534,28 @@ void PPPMDipole::compute(int eflag, int vflag) void PPPMDipole::allocate() { + // create ghost grid object for rho and electric field communication + // returns local owned and ghost grid bounds + // setup communication patterns and buffers + + gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm); + gc_dipole->set_distance(0.5*neighbor->skin + qdist); + gc_dipole->set_stencil_atom(-nlower,nupper); + gc_dipole->set_shift_atom(shiftatom); + gc_dipole->set_zfactor(slab_volfactor); + + gc_dipole->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); + + gc_dipole->setup_comm(ngc_buf1,ngc_buf2); + + npergrid = 9; + + memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1"); + memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2"); + + // allocate distributed grid data + memory->create3d_offset(densityx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm_dipole:densityx_brick_dipole"); memory->create3d_offset(densityy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, @@ -601,20 +629,6 @@ void PPPMDipole::allocate() nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, 1,0,0,FFT_PRECISION,collective_flag); - - // create ghost grid object for rho and electric field communication - // also create 2 bufs for ghost grid cell comm, passed to Grid3d methods - - gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); - - gc_dipole->setup(ngc_buf1,ngc_buf2); - - npergrid = 9; - - memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1"); - memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2"); } /* ---------------------------------------------------------------------- @@ -623,6 +637,10 @@ void PPPMDipole::allocate() void PPPMDipole::deallocate() { + delete gc_dipole; + memory->destroy(gc_buf1); + memory->destroy(gc_buf2); + memory->destroy3d_offset(densityx_brick_dipole,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(densityy_brick_dipole,nzlo_out,nylo_out,nxlo_out); memory->destroy3d_offset(densityz_brick_dipole,nzlo_out,nylo_out,nxlo_out); @@ -645,7 +663,9 @@ void PPPMDipole::deallocate() memory->destroy(work3); memory->destroy(work4); - delete gc_dipole; + delete fft1; + delete fft2; + delete remap; } /* ---------------------------------------------------------------------- diff --git a/src/KSPACE/pppm_dipole_spin.cpp b/src/KSPACE/pppm_dipole_spin.cpp index ea92eb4685..4352ac03ae 100644 --- a/src/KSPACE/pppm_dipole_spin.cpp +++ b/src/KSPACE/pppm_dipole_spin.cpp @@ -26,6 +26,7 @@ #include "grid3d.h" #include "math_const.h" #include "memory.h" +#include "neighbor.h" #include "pair.h" #include "update.h" @@ -173,7 +174,7 @@ void PPPMDipoleSpin::init() // or overlap is allowed, then done // else reduce order and try again - Grid3d *gctmp = nullptr; + gc_dipole = nullptr; int iteration = 0; while (order >= minorder) { @@ -186,23 +187,28 @@ void PPPMDipoleSpin::init() set_grid_local(); if (overlap_allowed) break; - gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, - nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, - nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); + gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm); + gc_dipole->set_distance(0.5*neighbor->skin + qdist); + gc_dipole->set_stencil_atom(-nlower,nupper); + gc_dipole->set_shift_atom(shiftatom); + gc_dipole->set_zfactor(slab_volfactor); + + gc_dipole->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); int tmp1,tmp2; - gctmp->setup(tmp1,tmp2); - if (gctmp->ghost_adjacent()) break; - delete gctmp; + gc_dipole->setup_comm(tmp1,tmp2); + if (gc_dipole->ghost_adjacent()) break; + delete gc_dipole; order--; iteration++; } if (order < minorder) error->all(FLERR,"PPPMDipoleSpin order < minimum allowed order"); - if (!overlap_allowed && !gctmp->ghost_adjacent()) + if (!overlap_allowed && !gc_dipole->ghost_adjacent()) error->all(FLERR,"PPPMDipoleSpin grid stencil extends beyond nearest neighbor processor"); - if (gctmp) delete gctmp; + if (gc_dipole) delete gc_dipole; // adjust g_ewald diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 1f73af81cc..336e7ca1c1 100644 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -446,7 +446,7 @@ void PPPMDisp::init() nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); int tmp1,tmp2; - gctmp->setup(tmp1,tmp2); + gctmp->setup_comm(tmp1,tmp2); if (gctmp->ghost_adjacent()) break; delete gctmp; @@ -526,7 +526,7 @@ void PPPMDisp::init() nzlo_out_6,nzhi_out_6); int tmp1,tmp2; - gctmp->setup(tmp1,tmp2); + gctmp->setup_comm(tmp1,tmp2); if (gctmp->ghost_adjacent()) break; delete gctmp; @@ -1754,7 +1754,7 @@ void _noopt PPPMDisp::allocate() nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out); - gc->setup(ngc_buf1,ngc_buf2); + gc->setup_comm(ngc_buf1,ngc_buf2); if (differentiation_flag) npergrid = 1; else npergrid = 3; @@ -1838,7 +1838,7 @@ void _noopt PPPMDisp::allocate() nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); - gc6->setup(ngc6_buf1,ngc6_buf2); + gc6->setup_comm(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 1; else npergrid6 = 3; @@ -2001,7 +2001,7 @@ void _noopt PPPMDisp::allocate() nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); - gc6->setup(ngc6_buf1,ngc6_buf2); + gc6->setup_comm(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 7; else npergrid6 = 21; @@ -2088,7 +2088,7 @@ void _noopt PPPMDisp::allocate() nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6, nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6); - gc6->setup(ngc6_buf1,ngc6_buf2); + gc6->setup_comm(ngc6_buf1,ngc6_buf2); if (differentiation_flag) npergrid6 = 1*nsplit_alloc; else npergrid6 = 3*nsplit_alloc; diff --git a/src/compute_property_grid.cpp b/src/compute_property_grid.cpp index 64092dd596..fec3920357 100644 --- a/src/compute_property_grid.cpp +++ b/src/compute_property_grid.cpp @@ -252,9 +252,9 @@ void ComputePropertyGrid::allocate_grid() ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1); } else { - grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, 0.0, 0, 0.0, nxlo_in, nxhi_in, nylo_in, - nyhi_in, nzlo_in, nzhi_in, nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, - nzhi_out); + grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid); + 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); if (nvalues == 1) memory->create3d_offset(vec3d, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out, "property/grid:vec3d"); diff --git a/src/dump_grid.cpp b/src/dump_grid.cpp index e8428ba5f2..72054145c1 100644 --- a/src/dump_grid.cpp +++ b/src/dump_grid.cpp @@ -509,7 +509,7 @@ int DumpGrid::count() grid3d = (Grid3d *) compute[field2index[0]]->get_grid_by_index(field2grid[0]); else if (field2source[0] == FIX) grid3d = (Grid3d *) fix[field2index[0]]->get_grid_by_index(field2grid[0]); - grid3d->get_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in); + grid3d->get_bounds_owned(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in); } // invoke Computes for per-grid quantities diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index 39b05fb105..508b0ddf51 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -411,10 +411,10 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1); } else { - grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 0, 0.5, - nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in, - nxlo_out, nxhi_out, nylo_out, nyhi_out, - nzlo_out, nzhi_out); + 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 diff --git a/src/grid3d.cpp b/src/grid3d.cpp index 56f52dfe5b..ac4cc6f91f 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -40,132 +40,29 @@ static constexpr int OFFSET = 16384; ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - constructor called by all classes except PPPM and MSM - comm_caller = caller's communicator - nx,ny,nz caller = size of global grid - maxdist = max distance outside of proc domain a particle will be - extra = additional ghost grid pts needed in each dim, e.g. for stencil - shift_caller = 0.0 for grid pt in lower-left corner of grid cell, - 0.5 for center - return: - i xyz lohi = portion of global grid this proc owns, 0 <= index < N - o xyz lohi = owned + ghost grid cells needed in all directions - for periodic dims, o indices can be < 0 or >= N - for non-periodic dims, o indices will be >= 0 and < N - since no grid comm is done across non-periodic boundaries - NOTE: allow zfactor to be a calling arg for PPPM ? + constructor to create a 3d distributed grid + gcomm = caller's communicator + gnx,gny,gnz = global grid size ------------------------------------------------------------------------- */ -Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm comm_caller, - int nx_caller, int ny_caller, int nz_caller, - double maxdist, int extra, double shift_caller, - int &ixlo, int &ixhi, int &iylo, int &iyhi, int &izlo, int &izhi, - int &oxlo, int &oxhi, int &oylo, int &oyhi, int &ozlo, int &ozhi) - : Pointers(lmp) +Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) : + Pointers(lmp) { - // store commnicator and global grid size - // set layout mode - - gridcomm = comm_caller; + gridcomm = gcomm; MPI_Comm_rank(gridcomm,&me); MPI_Comm_size(gridcomm,&nprocs); - nx = nx_caller; - ny = ny_caller; - nz = nz_caller; + nx = gnx; + ny = gny; + nz = gnz; + + // default settings, can be overridden by set() methods - // NOTE: hardwire shift = 0.5 ? - - shift = shift_caller; - - // define owned grid cells for each proc - // extend bounds with ghost grid cells in each direction - - double fraclo,frachi; - - if (comm->layout != Comm::LAYOUT_TILED) { - fraclo = comm->xsplit[comm->myloc[0]]; - frachi = comm->xsplit[comm->myloc[0]+1]; - partition_grid(nx,fraclo,frachi,shift,inxlo,inxhi); - fraclo = comm->ysplit[comm->myloc[1]]; - frachi = comm->ysplit[comm->myloc[1]+1]; - partition_grid(ny,fraclo,frachi,shift,inylo,inyhi); - fraclo = comm->zsplit[comm->myloc[2]]; - frachi = comm->zsplit[comm->myloc[2]+1]; - partition_grid(nz,fraclo,frachi,shift,inzlo,inzhi); - } else { - fraclo = comm->mysplit[0][0]; - frachi = comm->mysplit[0][1]; - partition_grid(nx,fraclo,frachi,shift,inxlo,inxhi); - fraclo = comm->mysplit[1][0]; - frachi = comm->mysplit[1][1]; - partition_grid(ny,fraclo,frachi,shift,inylo,inyhi); - fraclo = comm->mysplit[2][0]; - frachi = comm->mysplit[2][1]; - partition_grid(nz,fraclo,frachi,shift,inzlo,inzhi); - } - - ghost_grid(maxdist,extra); - - // error check on size of grid stored by this proc - - bigint total = (bigint) - (outxhi - outxlo + 1) * (outyhi - outylo + 1) * (outzhi - outzlo + 1); - if (total > MAXSMALLINT) - error->one(FLERR, "Too many owned+ghost grid3d points"); - - // default = caller grid is allocated to ghost grid - // used when computing pack/unpack lists in indices() - // NOTE: allow caller to override this - - fullxlo = outxlo; - fullxhi = outxhi; - fullylo = outylo; - fullyhi = outyhi; - fullzlo = outzlo; - fullzhi = outzhi; - - // initialize data structs - - nswap = maxswap = 0; - swap = nullptr; - - nsend = nrecv = ncopy = 0; - send = nullptr; - recv = nullptr; - copy = nullptr; - requests = nullptr; - requests_remap = nullptr; - - xsplit = ysplit = zsplit = nullptr; - grid2proc = nullptr; - rcbinfo = nullptr; - - nsend_remap = nrecv_remap = self_remap = 0; - send_remap = nullptr; - recv_remap = nullptr; - - // store info about Comm decomposition needed for remap operation - // two Grid instances will exist for duration of remap - // each must know Comm decomp at time Grid instance was created - - extract_comm_info(); - - // return values - - ixlo = inxlo; - ixhi = inxhi; - iylo = inylo; - iyhi = inyhi; - izlo = inzlo; - izhi = inzhi; - - oxlo = outxlo; - oxhi = outxhi; - oylo = outylo; - oyhi = outyhi; - ozlo = outzlo; - ozhi = outzhi; + maxdist = 0.0; + stencil_grid_lo = stencil_grid_hi = 0; + stencil_atom_lo = stencil_atom_hi = 0; + shift_grid = shift_atom = 0.0; + zfactor = 1.0; } /* ---------------------------------------------------------------------- @@ -367,7 +264,78 @@ void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi, } // ---------------------------------------------------------------------- -// access Grid parameters +// set Grid parameters +// ---------------------------------------------------------------------- + +/* ---------------------------------------------------------------------- + maxdist = max distance outside proc subdomain a particle can be +------------------------------------------------------------------------- */ + +void Grid3d::set_distance(double distance) +{ + maxdist = distance; +} + +/* ---------------------------------------------------------------------- + # of grid cells beyond an owned grid cell that caller accesses + e.g. for a finite different stencial + can be different in lo vs hi direction +------------------------------------------------------------------------- */ + +void Grid3d::set_stencil_grid(int lo, int hi) +{ + stencil_grid_lo = lo; + stencil_grid_hi = hi; +} + +/* ---------------------------------------------------------------------- + # of grid cells beyond a particle's grid cell that caller updates + e.g. for smearing a point charge to the grid + can be different in lo vs hi direction +------------------------------------------------------------------------- */ + +void Grid3d::set_stencil_atom(int lo, int hi) +{ + stencil_atom_lo = lo; + stencil_atom_hi = hi; +} + +/* ---------------------------------------------------------------------- + shift_grid = offset of position of grid point within grid cell + 0.5 = cell center, 0.0 = lower-left corner of cell + used to determine which proc owns the grid cell within its subdomain +------------------------------------------------------------------------- */ + +void Grid3d::set_shift_grid(double shift) +{ + shift_grid = shift; +} + +/* ---------------------------------------------------------------------- + shift_atom = offset of atoms when caller maps them to grid cells + 0.5 = half a grid cell, 0.0 = no offset + used when computing possible ghost extent + PPPM uses 0.5 when order is odd, 0.0 when order is even +------------------------------------------------------------------------- */ + +void Grid3d::set_shift_atom(double shift) +{ + shift_atom = shift; +} + +/* ---------------------------------------------------------------------- + zfactor > 1.0 when grid extends beyond simulation box bounds in Z + used by KSpace PPPM for 2d periodic slab geometries + default = 1.0 (no extension in Z) +------------------------------------------------------------------------- */ + +void Grid3d::set_zfactor(double factor) +{ + zfactor = factor; +} + +// ---------------------------------------------------------------------- +// retrieve Grid parameters // ---------------------------------------------------------------------- int Grid3d::identical(Grid3d *grid2) @@ -375,7 +343,7 @@ int Grid3d::identical(Grid3d *grid2) int inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2; int outxlo2,outxhi2,outylo2,outyhi2,outzlo2,outzhi2; - grid2->get_bounds(inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2); + grid2->get_bounds_owned(inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2); grid2->get_bounds_ghost(outxlo2,outxhi2,outylo2,outyhi2,outzlo2,outzhi2); int flag = 0; @@ -404,8 +372,8 @@ void Grid3d::get_size(int &nxgrid, int &nygrid, int &nzgrid) /* ---------------------------------------------------------------------- */ -void Grid3d::get_bounds(int &xlo, int &xhi, int &ylo, int &yhi, - int &zlo, int &zhi) +void Grid3d::get_bounds_owned(int &xlo, int &xhi, int &ylo, int &yhi, + int &zlo, int &zhi) { xlo = inxlo; xhi = inxhi; @@ -433,6 +401,113 @@ void Grid3d::get_bounds_ghost(int &xlo, int &xhi, int &ylo, int &yhi, // also store comm and grid partitioning info // ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- + setup grid partition for each proc = owned + ghost cells + return: + i xyz lohi = portion of global grid this proc owns, 0 <= index < N + o xyz lohi = owned + ghost grid cells in all directions + for periodic dims, o indices can be < 0 or >= N + for non-periodic dims, o indices will be >= 0 and < N + since no grid comm is done across non-periodic boundaries +------------------------------------------------------------------------- */ + +void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, + int &izlo, int &izhi, + int &oxlo, int &oxhi, int &oylo, int &oyhi, + 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) { + fraclo = comm->xsplit[comm->myloc[0]]; + frachi = comm->xsplit[comm->myloc[0]+1]; + partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi); + fraclo = comm->ysplit[comm->myloc[1]]; + frachi = comm->ysplit[comm->myloc[1]+1]; + partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + fraclo = comm->zsplit[comm->myloc[2]]; + frachi = comm->zsplit[comm->myloc[2]+1]; + partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi); + } else { + fraclo = comm->mysplit[0][0]; + frachi = comm->mysplit[0][1]; + partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi); + fraclo = comm->mysplit[1][0]; + frachi = comm->mysplit[1][1]; + partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + fraclo = comm->mysplit[2][0]; + frachi = comm->mysplit[2][1]; + partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi); + } + + // extend owned grid bounds with ghost grid cells in each direction + + ghost_grid(); + + // error check on size of grid stored by this proc + + bigint total = (bigint) + (outxhi - outxlo + 1) * (outyhi - outylo + 1) * (outzhi - outzlo + 1); + if (total > MAXSMALLINT) + error->one(FLERR, "Too many owned+ghost grid3d points"); + + // default = caller grid is allocated to ghost grid + // used when computing pack/unpack lists in indices() + // NOTE: allow caller to override this + + fullxlo = outxlo; + fullxhi = outxhi; + fullylo = outylo; + fullyhi = outyhi; + fullzlo = outzlo; + fullzhi = outzhi; + + // initialize data structs + + nswap = maxswap = 0; + swap = nullptr; + + nsend = nrecv = ncopy = 0; + send = nullptr; + recv = nullptr; + copy = nullptr; + requests = nullptr; + requests_remap = nullptr; + + xsplit = ysplit = zsplit = nullptr; + grid2proc = nullptr; + rcbinfo = nullptr; + + nsend_remap = nrecv_remap = self_remap = 0; + send_remap = nullptr; + recv_remap = nullptr; + + // store info about Comm decomposition needed for remap operation + // two Grid instances will exist for duration of remap + // each must know Comm decomp at time Grid instance was created + + extract_comm_info(); + + // return values + + ixlo = inxlo; + ixhi = inxhi; + iylo = inylo; + iyhi = inyhi; + izlo = inzlo; + izhi = inzhi; + + oxlo = outxlo; + oxhi = outxhi; + oylo = outylo; + oyhi = outyhi; + ozlo = outzlo; + ozhi = outzhi; +} + /* ---------------------------------------------------------------------- partition a global regular grid into one brick-shaped sub-grid per proc if grid point is inside my sub-domain I own it, @@ -468,7 +543,7 @@ void Grid3d::partition_grid(int ngrid, double fraclo, double frachi, ghost cell indices for periodic systems can be < 0 or >= N ------------------------------------------------------------------------- */ -void Grid3d::ghost_grid(double maxdist, int extra) +void Grid3d::ghost_grid() { double *prd,*boxlo,*sublo,*subhi; int triclinic = domain->triclinic; @@ -491,10 +566,11 @@ void Grid3d::ghost_grid(double maxdist, int extra) if (triclinic == 0) dist[0] = dist[1] = dist[2] = maxdist; else MathExtra::tribbox(domain->h,maxdist,&dist[0]); - // nlo,nhi = min/max index of global grid cell my owned atoms can be mapped to - // including up to maxdist displacement outside my subdomain - // extra = additional ghost layers required by called (e.g. finite diff stencil) - // max of the two quantities = ghost cell layers needed in each dim/dir + // lo/hi = min/max index of global grid cells my owned atoms can be mapped to + // includes effects of maxdist and shift_atom settings + // lo/hi can be further extended by stencil_atom and stencil_grid settings + // all those settings are set by caller + // ghost cell layers needed in each dim/dir = max of two extension effects // OFFSET allows generation of negative indices with static_cast // out xyz lo/hi = index range of owned + ghost cells @@ -503,22 +579,22 @@ void Grid3d::ghost_grid(double maxdist, int extra) double dzinv = nz / prd[2]; int lo, hi; - lo = static_cast((sublo[0]-dist[0]-boxlo[0]) * dxinv + OFFSET) - OFFSET; - hi = static_cast((subhi[0]+dist[0]-boxlo[0]) * dxinv + OFFSET) - OFFSET; - outxlo = MIN(lo, inxlo - extra); - outxhi = MAX(hi, inxhi + extra); + lo = static_cast((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom + OFFSET) - OFFSET; + hi = static_cast((subhi[0]+dist[0]-boxlo[0]) * dxinv + shift_atom + OFFSET) - OFFSET; + outxlo = MIN(lo - stencil_atom_lo, inxlo - stencil_grid_lo); + outxhi = MAX(hi + stencil_atom_hi, inxhi + stencil_grid_hi); - lo = static_cast((sublo[1]-dist[1]-boxlo[1]) * dyinv + OFFSET) - OFFSET; - hi = static_cast((subhi[1]+dist[1]-boxlo[1]) * dyinv + OFFSET) - OFFSET; - outylo = MIN(lo, inylo - extra); - outyhi = MAX(hi, inyhi + extra); + lo = static_cast((sublo[1]-dist[1]-boxlo[1]) * dyinv + shift_atom + OFFSET) - OFFSET; + hi = static_cast((subhi[1]+dist[1]-boxlo[1]) * dyinv + shift_atom + OFFSET) - OFFSET; + outxlo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo); + outxhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi); // NOTE: need to account for zfactor here - lo = static_cast((sublo[2]-dist[2]-boxlo[2]) * dzinv + OFFSET) - OFFSET; - hi = static_cast((subhi[2]+dist[2]-boxlo[2]) * dzinv + OFFSET) - OFFSET; - outzlo = MIN(lo, inzlo - extra); - outzhi = MAX(hi, inzhi + extra); + lo = static_cast((sublo[2]-dist[2]-boxlo[2]) * dzinv + shift_atom + OFFSET) - OFFSET; + hi = static_cast((subhi[2]+dist[2]-boxlo[2]) * dzinv + shift_atom + OFFSET) - OFFSET; + outxlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo); + outxhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi); // limit out xyz lo/hi indices to global grid for non-periodic dims @@ -1427,7 +1503,8 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2) // overlap_old = vector of overlap info in Overlap data struct int oldbox[6]; - old->get_bounds(oldbox[0],oldbox[1],oldbox[2],oldbox[3],oldbox[4],oldbox[5]); + old->get_bounds_owned(oldbox[0],oldbox[1],oldbox[2],oldbox[3], + oldbox[4],oldbox[5]); pbc[0] = pbc[1] = pbc[2] = 0; Overlap *overlap_old; @@ -1466,7 +1543,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2) // overlap_new = vector of overlap info in Overlap data struct int newbox[6]; - get_bounds(newbox[0],newbox[1],newbox[2],newbox[3],newbox[4],newbox[5]); + get_bounds_owned(newbox[0],newbox[1],newbox[2],newbox[3],newbox[4],newbox[5]); pbc[0] = pbc[1] = pbc[2] = 0; Overlap *overlap_new; @@ -1762,9 +1839,9 @@ int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap for (int k = kproclo; k <= kprochi; k++) for (int j = jproclo; j <= jprochi; j++) for (int i = iproclo; i <= iprochi; i++) { - partition_grid(nx,xsplit[i],xsplit[i+1],shift,obox[0],obox[1]); - partition_grid(ny,ysplit[j],ysplit[j+1],shift,obox[2],obox[3]); - partition_grid(nz,zsplit[k],zsplit[k+1],shift,obox[4],obox[5]); + partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,obox[0],obox[1]); + partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,obox[2],obox[3]); + partition_grid(nz,zsplit[k],zsplit[k+1],shift_grid,obox[4],obox[5]); if (noverlap_list == maxoverlap_list) grow_overlap(); overlap_list[noverlap_list].proc = grid2proc[i][j][k]; diff --git a/src/grid3d.h b/src/grid3d.h index c6047d2ddd..3eb3424a8b 100644 --- a/src/grid3d.h +++ b/src/grid3d.h @@ -22,21 +22,31 @@ class Grid3d : protected Pointers { public: enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes - Grid3d(class LAMMPS *, MPI_Comm, int, int, int, double, int, double, - int &, int &, int &, int &, int &, int &, - int &, int &, int &, int &, int &, int &); + Grid3d(class LAMMPS *, MPI_Comm, int, int, int); + Grid3d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int); Grid3d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int, int); + ~Grid3d() override; + void set_distance(double); + void set_stencil_grid(int, int); + void set_stencil_atom(int, int); + void set_shift_grid(double); + void set_shift_atom(double); + void set_zfactor(double); + int identical(Grid3d *); void get_size(int &, int &, int &); - void get_bounds(int &, int &, int &, int &, int &, int &); + void get_bounds_owned(int &, int &, int &, int &, int &, int &); void get_bounds_ghost(int &, int &, int &, int &, int &, int &); + void setup_grid(int &, int &, int &, int &, int &, int &, + int &, int &, int &, int &, int &, int &); + void setup_comm(int &, int &); int ghost_adjacent(); void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype); @@ -54,9 +64,21 @@ class Grid3d : protected Pointers { MPI_Comm gridcomm; // communicator for this class // usually world, but MSM calls with subset - // inputs from caller via constructor + // inputs from caller int nx, ny, nz; // size of global grid in all 3 dims + double maxdist; // distance owned atoms can move outside subdomain + int stencil_atom_lo,stencil_atom_hi; // grid cells accessed beyond atom's cell + int stencil_grid_lo,stencil_grid_hi; // grid cells accessed beyond owned cell + double shift_grid; // location of grid point within grid cell + // only affects which proc owns grid cell + double shift_atom; // shift applied to atomd when mapped to grid cell by caller + // only affects extent of ghost cells + double zfactor; // multiplier on extend of grid in Z direction + // 1.0 = no extra grid cells in Z + + // extent of my owned and ghost cells + int inxlo, inxhi; // inclusive extent of my grid chunk, 0 <= in <= N-1 int inylo, inyhi; int inzlo, inzhi; @@ -67,9 +89,6 @@ class Grid3d : protected Pointers { int fullylo, fullyhi; // can be same as out indices or larger int fullzlo, fullzhi; - double shift; // location of grid point within grid cell - // only affects which proc owns grid cell - // ------------------------------------------- // internal variables for BRICK layout // ------------------------------------------- @@ -215,7 +234,7 @@ class Grid3d : protected Pointers { // ------------------------------------------- void partition_grid(int, double, double, double, int &, int &); - void ghost_grid(double, int); + void ghost_grid(); void extract_comm_info(); void store(int, int, int, int, int, int, int, int, int, int, int, int,