From 1c5bdadcfb37dfbd05ef26afe24e1d03fd519934 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 18 Aug 2021 15:43:54 -0600 Subject: [PATCH] small alteration to code that assigns grid pts to procs --- src/KSPACE/msm.cpp | 19 +++----- src/KSPACE/pppm.cpp | 34 +++----------- src/KSPACE/pppm.h | 3 ++ src/KSPACE/pppm_disp.cpp | 18 ++------ src/comm.cpp | 97 ++++++++++++++++++++++++++++++++++++++++ src/comm.h | 6 +++ src/domain.h | 14 ++++-- 7 files changed, 133 insertions(+), 58 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 87b919af17..0a642d5b53 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -1194,18 +1194,13 @@ void MSM::set_grid_local() for (int n=0; n (comm->xsplit[comm->myloc[0]] * nx_msm[n]); - nxhi_in[n] = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_msm[n]) - 1; - - nylo_in[n] = static_cast (comm->ysplit[comm->myloc[1]] * ny_msm[n]); - nyhi_in[n] = static_cast (comm->ysplit[comm->myloc[1]+1] * ny_msm[n]) - 1; - - nzlo_in[n] = static_cast (comm->zsplit[comm->myloc[2]] * nz_msm[n]); - nzhi_in[n] = static_cast (comm->zsplit[comm->myloc[2]+1] * nz_msm[n]) - 1; + // partition global grid across procs + // nxyz lo/hi = lower/upper bounds of global grid this proc owns + // indices range from 0 to N-1 inclusive in each dim + + comm->partition_grid(nx_msm[n],ny_msm[n],nz_msm[n],0.0, + nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n], + nzlo_in[n],nzhi_in[n]); // nlower,nupper = stencil size for mapping (interpolating) particles to MSM grid diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 0415631a4c..1c77bd8919 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -1320,34 +1320,12 @@ double PPPM::final_accuracy() void PPPM::set_grid_local() { - // global indices of PPPM grid range from 0 to N-1 - // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of - // global PPPM grid that I own without ghost cells - // for slab PPPM, assign z grid as if it were not extended - // both non-tiled and tiled proc layouts use 0-1 fractional sumdomain info + // partition global grid across procs + // nxyz lo/hi = lower/upper bounds of global grid this proc owns + // indices range from 0 to N-1 inclusive in each dim - if (comm->layout != Comm::LAYOUT_TILED) { - nxlo_in = static_cast (comm->xsplit[comm->myloc[0]] * nx_pppm); - nxhi_in = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_pppm) - 1; - - nylo_in = static_cast (comm->ysplit[comm->myloc[1]] * ny_pppm); - nyhi_in = static_cast (comm->ysplit[comm->myloc[1]+1] * ny_pppm) - 1; - - nzlo_in = static_cast - (comm->zsplit[comm->myloc[2]] * nz_pppm/slab_volfactor); - nzhi_in = static_cast - (comm->zsplit[comm->myloc[2]+1] * nz_pppm/slab_volfactor) - 1; - - } else { - nxlo_in = static_cast (comm->mysplit[0][0] * nx_pppm); - nxhi_in = static_cast (comm->mysplit[0][1] * nx_pppm) - 1; - - nylo_in = static_cast (comm->mysplit[1][0] * ny_pppm); - nyhi_in = static_cast (comm->mysplit[1][1] * ny_pppm) - 1; - - nzlo_in = static_cast (comm->mysplit[2][0] * nz_pppm/slab_volfactor); - nzhi_in = static_cast (comm->mysplit[2][1] * nz_pppm/slab_volfactor) - 1; - } + comm->partition_grid(nx_pppm,ny_pppm,nz_pppm,slab_volfactor, + nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in); // nlower,nupper = stencil size for mapping particles to PPPM grid @@ -1367,7 +1345,7 @@ void PPPM::set_grid_local() // 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] = particle position bound = subbox + skin/2.0 + qdist + // 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 diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index 4fa41fa311..0b67040b00 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -120,6 +120,9 @@ class PPPM : public KSpace { double qdist; // distance from O site to negative charge double alpha; // geometric factor + void carve_grid(int, int, int, double, + int &, int &, int &, int &, int &, int &); + virtual void set_grid_global(); void set_grid_local(); void adjust_gewald(); diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index f7c97b1f87..a12c9ecd8f 100644 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -2700,22 +2700,10 @@ void PPPMDisp::set_fft_parameters(int& nx_p, int& ny_p, int& nz_p, int& ng, int& nf, int& nfb, double& sft, double& sftone, int& ord) { - // global indices of PPPM grid range from 0 to N-1 - // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of - // global PPPM grid that I own without ghost cells - // for slab PPPM, assign z grid as if it were not extended - - nxlo_i = static_cast (comm->xsplit[comm->myloc[0]] * nx_p); - nxhi_i = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_p) - 1; - - nylo_i = static_cast (comm->ysplit[comm->myloc[1]] * ny_p); - nyhi_i = static_cast (comm->ysplit[comm->myloc[1]+1] * ny_p) - 1; - - nzlo_i = static_cast - (comm->zsplit[comm->myloc[2]] * nz_p/slab_volfactor); - nzhi_i = static_cast - (comm->zsplit[comm->myloc[2]+1] * nz_p/slab_volfactor) - 1; + // partition global grid across procs + comm->partition_grid(nx_p,ny_p,nz_p,slab_volfactor, + nxlo_i,nxhi_i,nylo_i,nyhi_i,nzlo_i,nzhi_i); // nlow,nupp = stencil size for mapping particles to PPPM grid diff --git a/src/comm.cpp b/src/comm.cpp index 2dddba11c2..fa9175f62c 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -832,6 +832,103 @@ int Comm::binary(double value, int n, double *vec) return index; } +/* ---------------------------------------------------------------------- + partition a global regular grid into sub-grids matching proc sub-domains + nx,ny,nz = extent of global grid + indices into the global grid range from 0 to N-1 in each dim + zfactor = 0.0 if the grid exactly covers the simulation box + zfactor > 1.0 if the grid extends beyond the +z boundary by this factor + used by 2d slab-mode PPPM + this effectively maps proc sub-grids to a smaller subset of the grid + if grid point is inside my sub-domain I own it, + this includes sub-domain lo boundary but excludes hi boundary + nxyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + if proc owns no grid cells in a dim, then nlo > nhi + special case: 2 procs share boundary which a grid point is exactly on + 2 equality if tests insure a consistent decision as to which proc owns it +------------------------------------------------------------------------- */ + +void Comm::partition_grid(int nx, int ny, int nz, double zfactor, + int &nxlo, int &nxhi, int &nylo, int &nyhi, + int &nzlo, int &nzhi) +{ + double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi; + + if (layout != LAYOUT_TILED) { + xfraclo = xsplit[myloc[0]]; + xfrachi = xsplit[myloc[0]+1]; + yfraclo = ysplit[myloc[1]]; + yfrachi = ysplit[myloc[1]+1]; + zfraclo = zsplit[myloc[2]]; + zfrachi = zsplit[myloc[2]+1]; + } else { + xfraclo = mysplit[0][0]; + xfrachi = mysplit[0][1]; + yfraclo = mysplit[1][0]; + yfrachi = mysplit[1][1]; + zfraclo = mysplit[2][0]; + zfrachi = mysplit[2][1]; + } + + nxlo = static_cast (xfraclo * nx); + if (1.0*nxlo != xfraclo*nx) nxlo++; + nxhi = static_cast (xfrachi * nx); + if (1.0*nxhi == xfrachi*nx) nxhi--; + + nylo = static_cast (yfraclo * ny); + if (1.0*nylo != yfraclo*ny) nylo++; + nyhi = static_cast (yfrachi * ny); + if (1.0*nyhi == yfrachi*ny) nyhi--; + + if (zfactor == 0.0) { + nzlo = static_cast (zfraclo * nz); + if (1.0*nzlo != zfraclo*nz) nzlo++; + nzhi = static_cast (zfrachi * nz); + if (1.0*nzhi == zfrachi*nz) nzhi--; + } else { + nzlo = static_cast (zfraclo * nz/zfactor); + if (1.0*nzlo != zfraclo*nz) nzlo++; + nzhi = static_cast (zfrachi * nz/zfactor); + if (1.0*nzhi == zfrachi*nz) nzhi--; + } + + // OLD code + // could sometimes map grid points slightly outside a proc to the proc + + /* + if (layout != LAYOUT_TILED) { + nxlo = static_cast (xsplit[myloc[0]] * nx); + nxhi = static_cast (xsplit[myloc[0]+1] * nx) - 1; + + nylo = static_cast (ysplit[myloc[1]] * ny); + nyhi = static_cast (ysplit[myloc[1]+1] * ny) - 1; + + if (zfactor == 0.0) { + nzlo = static_cast (zsplit[myloc[2]] * nz); + nzhi = static_cast (zsplit[myloc[2]+1] * nz) - 1; + } else { + nzlo = static_cast (zsplit[myloc[2]] * nz/zfactor); + nzhi = static_cast (zsplit[myloc[2]+1] * nz/zfactor) - 1; + } + + } else { + nxlo = static_cast (mysplit[0][0] * nx); + nxhi = static_cast (mysplit[0][1] * nx) - 1; + + nylo = static_cast (mysplit[1][0] * ny); + nyhi = static_cast (mysplit[1][1] * ny) - 1; + + if (zfactor == 0.0) { + nzlo = static_cast (mysplit[2][0] * nz); + nzhi = static_cast (mysplit[2][1] * nz) - 1; + } else { + nzlo = static_cast (mysplit[2][0] * nz/zfactor); + nzhi = static_cast (mysplit[2][1] * nz/zfactor) - 1; + } + } + */ +} + /* ---------------------------------------------------------------------- communicate inbuf around full ring of processors with messtag nbytes = size of inbuf = n datums * nper bytes diff --git a/src/comm.h b/src/comm.h index bc5faa49f4..7ad9e8756c 100644 --- a/src/comm.h +++ b/src/comm.h @@ -105,6 +105,11 @@ class Comm : protected Pointers { virtual void coord2proc_setup() {} virtual int coord2proc(double *, int &, int &, int &); + // partition a global regular grid by proc sub-domains + + void partition_grid(int, int, int, double, + int &, int &, int &, int &, int &, int &); + // memory usage virtual double memory_usage() = 0; @@ -117,6 +122,7 @@ class Comm : protected Pointers { int statflag = 0); // extract data useful to other classes + virtual void *extract(const char *, int &) { return nullptr; } protected: diff --git a/src/domain.h b/src/domain.h index 2a812218f1..f05afc50e6 100644 --- a/src/domain.h +++ b/src/domain.h @@ -42,38 +42,46 @@ class Domain : protected Pointers { int tiltsmall; // 1 if limit tilt, else 0 // orthogonal box + double xprd, yprd, zprd; // global box dimensions double xprd_half, yprd_half, zprd_half; // half dimensions double prd[3]; // array form of dimensions double prd_half[3]; // array form of half dimensions // triclinic box - // xprd,xprd_half,prd,prd_half = - // same as if untilted + // xyzprd,xyzprd_half and prd,prd_half = same as if untilted + double prd_lamda[3]; // lamda box = (1,1,1) double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5) - double boxlo[3], boxhi[3]; // orthogonal box global bounds + // orthogonal box global bounds + + double boxlo[3], boxhi[3]; // triclinic box // boxlo/hi = same as if untilted + double boxlo_lamda[3], boxhi_lamda[3]; // lamda box = (0,1) double boxlo_bound[3], boxhi_bound[3]; // bounding box of tilted domain double corners[8][3]; // 8 corner points // orthogonal box & triclinic box + double minxlo, minxhi; // minimum size of global box double minylo, minyhi; // when shrink-wrapping double minzlo, minzhi; // tri only possible for non-skew dims // orthogonal box + double sublo[3], subhi[3]; // sub-box bounds on this proc // triclinic box // sublo/hi = undefined + double sublo_lamda[3], subhi_lamda[3]; // bounds of subbox in lamda // triclinic box + double xy, xz, yz; // 3 tilt factors double h[6], h_inv[6]; // shape matrix in Voigt ordering // Voigt = xx,yy,zz,yz,xz,xy