From 8af384243f84b0e8b0995cf709490edf4dd85cf3 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 2 Nov 2022 15:01:58 -0600 Subject: [PATCH] support for more caller options in Grid2d/3d --- src/KSPACE/msm.cpp | 81 ++++------- src/grid2d.cpp | 344 ++++++++++++++++++++++----------------------- src/grid2d.h | 21 +-- src/grid3d.cpp | 224 ++++++++++++++++++----------- src/grid3d.h | 12 +- 5 files changed, 357 insertions(+), 325 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 0705a91e64..60930943b0 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -594,32 +594,17 @@ void MSM::allocate() memory->create2d_offset(phi1d,3,-order,order,"msm:phi1d"); memory->create2d_offset(dphi1d,3,-order,order,"msm:dphi1d"); - // one Grid3d using all processors for finest grid level - - gcall = new Grid3d(lmp,world,nx_msm[0],ny_msm[0],nz_msm[0]); - gcall->set_distance(0.5*neighbor->skin); - gcall->set_stencil_atom(-nlower,nupper); - - gcall->setup_grid(nxlo_in[0],nxhi_in[0],nylo_in[0], - nyhi_in[0],nzlo_in[0],nzhi_in[0], - nxlo_out_all,nxhi_out_all,nylo_out_all, - nyhi_out_all,nzlo_out_all,nzhi_out_all); - - gcall->set_larger_grid(nxlo_out[0],nxhi_out[0],nylo_out[0], - nyhi_out[0],nzlo_out[0],nzhi_out[0]); + // one Grid3d for finest grid level, using world comm and all procs + // use set_caller_grid() b/c MSM allocates local grid > out_all values - // NOTE: or is it out[0] ?? - // NOTE: worry about flag = 1 here, 2 later - - /* - gcall = new Grid3d(lmp,world,1,nx_msm[0],ny_msm[0],nz_msm[0], + gcall = new Grid3d(lmp,world,nx_msm[0],ny_msm[0],nz_msm[0], nxlo_in[0],nxhi_in[0],nylo_in[0], nyhi_in[0],nzlo_in[0],nzhi_in[0], nxlo_out_all,nxhi_out_all,nylo_out_all, - nyhi_out_all,nzlo_out_all,nzhi_out_all, - nxlo_out[0],nxhi_out[0],nylo_out[0], - nyhi_out[0],nzlo_out[0],nzhi_out[0]); - */ + nyhi_out_all,nzlo_out_all,nzhi_out_all); + + gcall->set_caller_grid(nxlo_out[0],nxhi_out[0],nylo_out[0], + nyhi_out[0],nzlo_out[0],nzhi_out[0]); gcall->setup_comm(ngcall_buf1,ngcall_buf2); npergrid = 1; @@ -639,34 +624,20 @@ void MSM::allocate() memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n], nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid"); - // create commgrid object for rho and electric field communication + // one Grid3d per level, using level-specific comm for coarser grids + // use set_proc_neigh() b/c MSM excludes non-participating procs from owned/ghost comm if (active_flag[n]) { delete gc[n]; - // NOTE: why is n = 0 same as all for grid size ? - - gc[n] = new Grid3d(lmp,world_levels[n],nx_msm[n],ny_msm[n],nz_msm[n]); - gc[n]->set_stencil_atom(-nlower,nupper); - - gc[n]->setup_grid(nxlo_in[n],nxhi_in[n],nylo_in[n], - nyhi_in[n],nzlo_in[n],nzhi_in[n], - nxlo_out[n],nxhi_out[n],nylo_out[n], - nyhi_out[n],nzlo_out[n],nzhi_out[n]); + gc[n] = new Grid3d(lmp,world_levels[n],nx_msm[n],ny_msm[n],nz_msm[n], + nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],nzlo_in[n],nzhi_in[n], + nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n],nzlo_out[n],nzhi_out[n]); int **procneigh = procneigh_levels[n]; gc[n]->set_proc_neighs(procneigh[0][0],procneigh[0][1],procneigh[1][0], procneigh[1][1],procneigh[2][0],procneigh[2][1]); - /* - gc[n] = new Grid3d(lmp,world_levels[n],2,nx_msm[n],ny_msm[n],nz_msm[n], - nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n], - nzlo_in[n],nzhi_in[n], - nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n], - procneigh[0][0],procneigh[0][1],procneigh[1][0], - procneigh[1][1],procneigh[2][0],procneigh[2][1]); - */ - gc[n]->setup_comm(ngc_buf1[n],ngc_buf2[n]); npergrid = 1; memory->destroy(gc_buf1[n]); @@ -1194,18 +1165,23 @@ void MSM::set_grid_local() for (int n = 0; n < levels; n++) { - // deleted and nullify grid arrays since the number or offset of gridpoints may change + // delete and nullify grid arrays since the number or offset of gridpoints may change memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); - // 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 + // use Grid3d to partition each grid into owned cells on each proc + // assignment of ghost cells is done below + // Grid3d is called later in allocate() with owned+ghost bounds - //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]); + gcall = new Grid3d(lmp,world,nx_msm[n],ny_msm[n],nz_msm[n]); + gcall->setup_grid(nxlo_in[n],nxhi_in[n],nylo_in[n], + nyhi_in[n],nzlo_in[n],nzhi_in[n], + nxlo_out[n],nxhi_out[n],nylo_out[n], + nyhi_out[n],nzlo_out[n],nzhi_out[n]); + delete gcall; + + // nlower/nupper = stencil size for mapping particles to grid nlower = -(order-1)/2; nupper = order/2; @@ -1236,10 +1212,12 @@ void MSM::set_grid_local() // position a particle in my box can be at // dist[3] = particle position bound = subbox + skin/2.0 // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping + // for n = 0, use a smaller ghost region for interpolation + // for n > 0, larger ghost region needed for direct sum and restriction/prolongation double dist[3]; double cuthalf = 0.0; - if (n == 0) cuthalf = 0.5*neighbor->skin; // only applies to finest grid + if (n == 0) cuthalf = 0.5*neighbor->skin; // only applies to finest grid dist[0] = dist[1] = dist[2] = cuthalf; if (triclinic) MathExtra::tribbox(domain->h,cuthalf,&dist[0]); @@ -1249,14 +1227,12 @@ void MSM::set_grid_local() nx_msm[n]/xprd + OFFSET) - OFFSET; nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * nx_msm[n]/xprd + OFFSET) - OFFSET; + if (n == 0) { - // use a smaller ghost region for interpolation nxlo_out_all = nlo + nlower; nxhi_out_all = nhi + nupper; } - // a larger ghost region is needed for the direct sum and for restriction/prolongation - nxlo_out[n] = nlo + MIN(-order,nxlo_direct); nxhi_out[n] = nhi + MAX(order,nxhi_direct); @@ -1268,6 +1244,7 @@ void MSM::set_grid_local() nylo_out_all = nlo + nlower; nyhi_out_all = nhi + nupper; } + nylo_out[n] = nlo + MIN(-order,nylo_direct); nyhi_out[n] = nhi + MAX(order,nyhi_direct); diff --git a/src/grid2d.cpp b/src/grid2d.cpp index 6a9385d30b..77ae289d63 100644 --- a/src/grid2d.cpp +++ b/src/grid2d.cpp @@ -40,9 +40,11 @@ static constexpr int OFFSET = 16384; ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - constructor to create a 3d distributed grid + constructor to create a 2d distributed grid + Grid3d assigns owned/ghost cells to each proc via setup_grid() + it MUST be called after constructor gcomm = caller's communicator - gnx,gny = global grid size + gnx,gny,gnz = global grid size ------------------------------------------------------------------------- */ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) : Pointers(lmp) @@ -55,32 +57,35 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) : Pointers(lmp) ny = gny; // default settings, can be overridden by set() methods + // these affect assignment of owned and ghost cells maxdist = 0.0; stencil_grid_lo = stencil_grid_hi = 0; stencil_atom_lo = stencil_atom_hi = 0; - shift_grid = shift_atom = 0.0; + shift_grid = 0.5; + shift_atom_lo = shift_atom_hi = 0.0; + yextra = 0; + yfactor = 1.0; } /* ---------------------------------------------------------------------- - constructor called by PPPM classes - gcomm = world communicator - gnx, gny = size of global grid - i xy lohi = portion of global grid this proc owns, 0 <= index < N - o xy lohi = owned grid portion + ghost grid cells needed in all directions - if o indices are < 0 or hi indices are >= N, - then grid is treated as periodic in that dimension, - communication is done across the periodic boundaries + alternate constructor to create a 2d distributed grid + caller assigns owned/ghost cells to each proc + setup_grid() must NOT be called + gcomm = caller's communicator + gnx,gny = global grid size + i xy lo/hi = extent of owned grid cells on this proc + o xy lo/hi = extent of owned+ghost grid cells on this proc + owned and ghost indices are inclusive + owned indices range from 0 to N-1 + ghost indices can extend < 0 or >= N ------------------------------------------------------------------------- */ -Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, - int gnx, int gny, +Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int ixlo, int ixhi, int iylo, int iyhi, - int oxlo, int oxhi, int oylo, int oyhi) - : Pointers(lmp) + int oxlo, int oxhi, int oylo, int oyhi) : + Pointers(lmp) { - // store commnicator and global grid size - gridcomm = gcomm; MPI_Comm_rank(gridcomm,&me); MPI_Comm_size(gridcomm,&nprocs); @@ -88,87 +93,22 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, nx = gnx; ny = gny; - // 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 + // store owned/ghost indices provided by caller + + inxlo = ixlo; + inxhi = ixhi; + inylo = iylo; + inyhi = iyhi; - extract_comm_info(); + outxlo = oxlo; + outxhi = oxhi; + outylo = oylo; + outyhi = oyhi; - // store grid bounds and proc neighs + // additional intialization + // other constructor invokes this from setup_grid() - if (layout != Comm::LAYOUT_TILED) { - store(ixlo,ixhi,iylo,iyhi, - oxlo,oxhi,oylo,oyhi, - oxlo,oxhi,oylo,oyhi, - procxlo,procxhi,procylo,procyhi); - } else { - store(ixlo,ixhi,iylo,iyhi, - oxlo,oxhi,oylo,oyhi, - oxlo,oxhi,oylo,oyhi, - 0,0,0,0); - } -} - -/* ---------------------------------------------------------------------- - constructor called by MSM - gcomm = world communicator or sub-communicator for a hierarchical grid - flag = 1 if e xy lohi values = larger grid stored by caller in gcomm = world - flag = 2 if e xy lohi values = 6 neighbor procs in gcomm - gnx, gny = size of global grid - i xy lohi = portion of global grid this proc owns, 0 <= index < N - o xy lohi = owned grid portion + ghost grid cells needed in all directions - e xy lohi for flag = 1: extent of larger grid stored by caller - e xy lohi for flag = 2: 4 neighbor procs -------------------------------------------------------------------------- */ - -Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag, - int gnx, int gny, - int ixlo, int ixhi, int iylo, int iyhi, - int oxlo, int oxhi, int oylo, int oyhi, - int exlo, int exhi, int eylo, int eyhi) - : Pointers(lmp) -{ - // store commnicator and global grid size - - gridcomm = gcomm; - MPI_Comm_rank(gridcomm,&me); - MPI_Comm_size(gridcomm,&nprocs); - - nx = gnx; - ny = gny; - - // 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(); - - // store grid bounds and proc neighs - - if (flag == 1) { - if (layout != Comm::LAYOUT_TILED) { - // this assumes gcomm = world - store(ixlo,ixhi,iylo,iyhi, - oxlo,oxhi,oylo,oyhi, - exlo,exhi,eylo,eyhi, - procxlo,procxhi,procylo,procyhi); - } else { - store(ixlo,ixhi,iylo,iyhi, - oxlo,oxhi,oylo,oyhi, - exlo,exhi,eylo,eyhi, - 0,0,0,0); - } - - } else if (flag == 2) { - if (layout != Comm::LAYOUT_TILED) { - store(ixlo,ixhi,iylo,iyhi, - oxlo,oxhi,oylo,oyhi, - oxlo,oxhi,oylo,oyhi, - exlo,exhi,eylo,eyhi); - } else { - error->all(FLERR,"Grid2d does not support tiled layout for MSM"); - } - } + initialize(); } /* ---------------------------------------------------------------------- */ @@ -213,42 +153,13 @@ Grid2d::~Grid2d() deallocate_remap(); } -/* ---------------------------------------------------------------------- - store grid bounds and proc neighs from caller in internal variables -------------------------------------------------------------------------- */ - -void Grid2d::store(int ixlo, int ixhi, int iylo, int iyhi, - int oxlo, int oxhi, int oylo, int oyhi, - int fxlo, int fxhi, int fylo, int fyhi, - int pxlo, int pxhi, int pylo, int pyhi) -{ - inxlo = ixlo; - inxhi = ixhi; - inylo = iylo; - inyhi = iyhi; - - outxlo = oxlo; - outxhi = oxhi; - outylo = oylo; - outyhi = oyhi; - - fullxlo = fxlo; - fullxhi = fxhi; - fullylo = fylo; - fullyhi = fyhi; - - procxlo = pxlo; - procxhi = pxhi; - procylo = pylo; - procyhi = pyhi; -} - // ---------------------------------------------------------------------- // set Grid parameters // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- maxdist = max distance outside proc subdomain a particle can be + used to determine extent of ghost cells ------------------------------------------------------------------------- */ void Grid2d::set_distance(double distance) @@ -258,7 +169,6 @@ void Grid2d::set_distance(double 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 ------------------------------------------------------------------------- */ @@ -269,8 +179,7 @@ void Grid2d::set_stencil_grid(int lo, int hi) } /* ---------------------------------------------------------------------- - # of grid cells beyond a particle's grid cell that caller updates - e.g. for smearing a point charge to the grid + # of grid cells beyond a particle's grid cell that caller accesses can be different in lo vs hi direction ------------------------------------------------------------------------- */ @@ -281,9 +190,9 @@ void Grid2d::set_stencil_atom(int lo, int 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 + shift_grid = offset within grid cell of position of grid point + 0.5 = cell center (default), 0.0 = lower-left corner of cell + used to decide which proc owns a grid cell (grid pt within subdomain) ------------------------------------------------------------------------- */ void Grid2d::set_shift_grid(double shift) @@ -292,15 +201,58 @@ void Grid2d::set_shift_grid(double shift) } /* ---------------------------------------------------------------------- - shift_atom = offset of atoms when caller maps them to grid cells + shift_atom = offset added to 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 + used to compute maximum possible ghost extents + use of lo/hi allows max ghost extent on each side to be different ------------------------------------------------------------------------- */ -void Grid2d::set_shift_atom(double shift) +void Grid2d::set_shift_atom(double shift_lo, double shift_hi) { - shift_atom = shift; + shift_atom_lo = shift_lo; + shift_atom_hi = shift_hi; +} + +/* ---------------------------------------------------------------------- + enable extra grid cells in Y + factor = muliplication factor on box size Y and thus grid size + factor > 1.0 when grid extends beyond Y box size (3.0 = tripled in size) + default zextra = 0, factor = 1.0 (no extra grid cells in Y) +------------------------------------------------------------------------- */ + +void Grid2d::set_yfactor(double factor) +{ + if (factor == 1.0) yextra = 0; + else yextra = 1; + yfactor = factor; +} + +/* ---------------------------------------------------------------------- + set IDs of proc neighbors used in uniform local owned/ghost comm + called AFTER setup_grid() but BEFORE setup_comm() to override + the processor neighbors stored by extract_comm() +------------------------------------------------------------------------- */ + +void Grid2d::set_proc_neighs(int pxlo, int pxhi, int pylo, int pyhi) +{ + procxlo = pxlo; + procxhi = pxhi; + procylo = pylo; + procyhi = pyhi; +} + +/* ---------------------------------------------------------------------- + set allocation dimensions of caller grid used by indices() to setup pack/unpack + called AFTER setup_grid() but BEFORE setup_comm() to override + the caller grid size set by setup_grid() and used in indices() +------------------------------------------------------------------------- */ + +void Grid2d::set_caller_grid(int fxlo, int fxhi, int fylo, int fyhi) +{ + fullxlo = fxlo; + fullxhi = fxhi; + fullylo = fylo; + fullyhi = fyhi; } // ---------------------------------------------------------------------- @@ -382,32 +334,56 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, 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); + partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi); fraclo = comm->ysplit[comm->myloc[1]]; frachi = comm->ysplit[comm->myloc[1]+1]; - partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + partition_grid(ny,fraclo,frachi,shift_grid,yextra,inylo,inyhi); } else { fraclo = comm->mysplit[0][0]; frachi = comm->mysplit[0][1]; - partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi); + partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi); fraclo = comm->mysplit[1][0]; frachi = comm->mysplit[1][1]; - partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + partition_grid(ny,fraclo,frachi,shift_grid,yextra,inylo,inyhi); } // extend owned grid bounds with ghost grid cells in each direction ghost_grid(); + // additional intialization + // other constructor invokes this directly + + initialize(); + + // return values + + ixlo = inxlo; + ixhi = inxhi; + iylo = inylo; + iyhi = inyhi; + + oxlo = outxlo; + oxhi = outxhi; + oylo = outylo; + oyhi = outyhi; +} + +/* ---------------------------------------------------------------------- + additional one-time setup common to both constructors + ---------------------------------------------------------------------- */ + +void Grid2d::initialize() +{ // error check on size of grid stored by this proc bigint total = (bigint) (outxhi - outxlo + 1) * (outyhi - outylo + 1); if (total > MAXSMALLINT) - error->one(FLERR, "Too many owned+ghost grid3d points"); + error->one(FLERR, "Too many owned+ghost grid2d points"); // default = caller grid is allocated to ghost grid // used when computing pack/unpack lists in indices() - // NOTE: allow caller to override this + // these values can be overridden using set_caller_grid() fullxlo = outxlo; fullxhi = outxhi; @@ -439,18 +415,6 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, // each must know Comm decomp at time Grid instance was created extract_comm_info(); - - // return values - - ixlo = inxlo; - ixhi = inxhi; - iylo = inylo; - iyhi = inyhi; - - oxlo = outxlo; - oxhi = outxhi; - oylo = outylo; - oyhi = outyhi; } /* ---------------------------------------------------------------------- @@ -460,11 +424,10 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, 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 - // NOTE: for future support of zfactor - 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 + 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 @@ -474,12 +437,19 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, ------------------------------------------------------------------------- */ void Grid2d::partition_grid(int ngrid, double fraclo, double frachi, - double shift, int &lo, int &hi) + double shift, int extra, int &lo, int &hi) { - lo = static_cast (fraclo * ngrid); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; - hi = static_cast (frachi * ngrid); - while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + if (extra == 0) { + lo = static_cast (fraclo * ngrid); + while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + hi = static_cast (frachi * ngrid); + while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + } else { + lo = static_cast (fraclo * ngrid/yfactor); + while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + hi = static_cast (frachi * ngrid/yfactor); + while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + } } /* ---------------------------------------------------------------------- @@ -521,19 +491,47 @@ void Grid2d::ghost_grid() double dxinv = nx / prd[0]; double dyinv = ny / prd[1]; + double dyinv_slab = ny / (prd[1] * yfactor); + int lo, hi; - 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; + lo = static_cast((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom_lo + OFFSET) - OFFSET; + hi = static_cast((subhi[0]+dist[0]-boxlo[0]) * dxinv + shift_atom_hi + 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 + shift_atom + OFFSET) - OFFSET; - hi = static_cast((subhi[1]+dist[1]-boxlo[1]) * dyinv + shift_atom + OFFSET) - OFFSET; - outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo); - outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi); + if (yextra == 0) { + lo = static_cast((sublo[1]-dist[1]-boxlo[1]) * dyinv + shift_atom_lo + OFFSET) - OFFSET; + hi = static_cast((subhi[1]+dist[1]-boxlo[1]) * dyinv + shift_atom_hi + OFFSET) - OFFSET; + outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo); + outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi); + } else { + lo = static_cast((sublo[1]-dist[1]-boxlo[1]) * dyinv_slab + shift_atom_lo + OFFSET) - OFFSET; + hi = static_cast((subhi[1]+dist[1]-boxlo[1]) * dyinv_slab + shift_atom_hi + OFFSET) - OFFSET; + outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo); + outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi); + } + + // if yextra = 1: + // adjust grid boundaries for processors at +y end, + // to include added empty grid cells between periodically repeating slabs + // in this case: + // want grid data forward communicated from +y proc to -y proc, but not vice versa + // want grid data reverse communicated from -y proc to +y proc, but not vice versa + // this is accomplished by inyhi = outyhi on +y end (no ghost cells) + // also insure no other procs use ghost cells beyond +y limit + + if (yextra) { + if (layout != Comm::LAYOUT_TILED) { + if (comm->myloc[1] == comm->procgrid[1]-1) inyhi = outyhi = ny - 1; + } else { + if (comm->mysplit[1][1] == 1.0) inyhi = outyhi = ny - 1; + } + outyhi = MIN(outyhi,ny-1); + } // limit out xyz lo/hi indices to global grid for non-periodic dims + // if yextra = 1, grid is still fully periodic int *periodicity = domain->periodicity; @@ -558,7 +556,7 @@ void Grid2d::extract_comm_info() layout = comm->layout; // for non TILED layout: - // proc xyz lohi = my 6 neighbor procs in this MPI_Comm + // proc xyz lohi = my 64neighbor procs in this MPI_Comm // NOTE: will need special logic for MSM case with different MPI_Comm // xyz split = copy of 1d vectors in Comm // grid2proc = copy of 3d array in Comm @@ -1651,8 +1649,8 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap for (int j = jproclo; j <= jprochi; j++) for (int i = iproclo; i <= iprochi; i++) { - 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(nx,xsplit[i],xsplit[i+1],shift_grid,0,obox[0],obox[1]); + partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,yextra,obox[2],obox[3]); if (noverlap_list == maxoverlap_list) grow_overlap(); overlap_list[noverlap_list].proc = grid2proc[i][j][0]; diff --git a/src/grid2d.h b/src/grid2d.h index d8e804c945..3ad28b3a4f 100644 --- a/src/grid2d.h +++ b/src/grid2d.h @@ -23,18 +23,17 @@ class Grid2d : protected Pointers { enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes Grid2d(class LAMMPS *, MPI_Comm, int, int); - Grid2d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int); - Grid2d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int, - int, int, int, int, int); - ~Grid2d() 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_shift_atom(double, double); + void set_yfactor(double); + void set_caller_grid(int, int, int, int); + void set_proc_neighs(int, int, int, int); int identical(Grid2d *); void get_size(int &, int &); @@ -68,8 +67,12 @@ protected: 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 + double shift_atom_lo,shift_atom_hi;; // max shift applied to atoms + // when mapped to grid cell by caller + // can be different in lo/hi directions // only affects extent of ghost cells + int yextra; // 1 if extra grid cells in Y, 0 if not + double yfactor; // multiplier on extent of grid in Y direction // extent of my owned and ghost cells @@ -224,13 +227,11 @@ protected: // internal methods // ------------------------------------------- - void partition_grid(int, double, double, double, int &, int &); + void initialize(); + void partition_grid(int, double, double, double, int, int &, int &); void ghost_grid(); void extract_comm_info(); - void store(int, int, int, int, int, int, int, int, - int, int, int, int, int, int, int, int); - void setup_comm_brick(int &, int &); void setup_comm_tiled(int &, int &); int ghost_adjacent_brick(); diff --git a/src/grid3d.cpp b/src/grid3d.cpp index 6dd4bd0110..30056d27c9 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -41,6 +41,8 @@ static constexpr int OFFSET = 16384; /* ---------------------------------------------------------------------- constructor to create a 3d distributed grid + Grid3d assigns owned/ghost cells to each proc via setup_grid() + it MUST be called after constructor gcomm = caller's communicator gnx,gny,gnz = global grid size ------------------------------------------------------------------------- */ @@ -57,6 +59,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) : nz = gnz; // default settings, can be overridden by set() methods + // these affect assignment of owned and ghost cells maxdist = 0.0; stencil_grid_lo = stencil_grid_hi = 0; @@ -64,6 +67,56 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) : shift_grid = 0.5; shift_atom_lo = shift_atom_hi = 0.0; zextra = 0; + zfactor = 1.0; +} + +/* ---------------------------------------------------------------------- + alternate constructor to create a 3d distributed grid + caller assigns owned/ghost cells to each proc + setup_grid() must NOT be called + used by MSM b/c its definition of ghost cells is complex + gcomm = caller's communicator + gnx,gny,gnz = global grid size + i xyz lo/hi = extent of owned grid cells on this proc + o xyz lo/hi = extent of owned+ghost grid cells on this proc + owned and ghost indices are inclusive + owned indices range from 0 to N-1 + ghost indices can extend < 0 or >= N +------------------------------------------------------------------------- */ + +Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz, + 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) +{ + gridcomm = gcomm; + MPI_Comm_rank(gridcomm,&me); + MPI_Comm_size(gridcomm,&nprocs); + + nx = gnx; + ny = gny; + nz = gnz; + + // store owned/ghost indices provided by caller + + inxlo = ixlo; + inxhi = ixhi; + inylo = iylo; + inyhi = iyhi; + inzlo = izlo; + inzhi = izhi; + + outxlo = oxlo; + outxhi = oxhi; + outylo = oylo; + outyhi = oyhi; + outzlo = ozlo; + outzhi = ozhi; + + // additional intialization + // other constructor invokes this from setup_grid() + + initialize(); } /* ---------------------------------------------------------------------- */ @@ -150,7 +203,7 @@ void Grid3d::set_stencil_atom(int lo, int hi) /* ---------------------------------------------------------------------- shift_grid = offset within grid cell of position of grid point 0.5 = cell center (default), 0.0 = lower-left corner of cell - used by Grid to decide which proc owns a grid cell (grid pt within subdomain) + used to decide which proc owns a grid cell (grid pt within subdomain) ------------------------------------------------------------------------- */ void Grid3d::set_shift_grid(double shift) @@ -161,7 +214,7 @@ void Grid3d::set_shift_grid(double shift) /* ---------------------------------------------------------------------- shift_atom = offset added to atoms when caller maps them to grid cells 0.5 = half a grid cell, 0.0 = no offset - used by Grid to compute maximum possible ghost extents + used to compute maximum possible ghost extents use of lo/hi allows max ghost extent on each side to be different PPPM uses 0.5 when stencil order is odd, 0.0 when order is even PPPM/stagger applies different shift values for 2 stagger iterations @@ -190,7 +243,7 @@ void Grid3d::set_zfactor(double factor) /* ---------------------------------------------------------------------- set IDs of proc neighbors used in uniform local owned/ghost comm - call this AFTER setup_grid() but BEFORE setup_comm() to override + called AFTER setup_grid() but BEFORE setup_comm() to override the processor neighbors stored by extract_comm() used by MSM to exclude non-participating procs for coarse grid comm ------------------------------------------------------------------------- */ @@ -208,13 +261,13 @@ void Grid3d::set_proc_neighs(int pxlo, int pxhi, int pylo, int pyhi, /* ---------------------------------------------------------------------- set allocation dimensions of caller grid used by indices() to setup pack/unpack - call this AFTER setup_grid() but BEFORE setup_comm() to override - the caller grid size set by setup_grid() + called AFTER setup_grid() but BEFORE setup_comm() to override + the caller grid size set by setup_grid() and used in indices() used by MSM to allow a larger level 0 grid to be allocated with more ghost cells for other operations ------------------------------------------------------------------------- */ -void Grid3d::set_larger_grid(int fxlo, int fxhi, int fylo, int fyhi, +void Grid3d::set_caller_grid(int fxlo, int fxhi, int fylo, int fyhi, int fzlo, int fzhi) { fullxlo = fxlo; @@ -308,36 +361,63 @@ 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) { fraclo = comm->xsplit[comm->myloc[0]]; frachi = comm->xsplit[comm->myloc[0]+1]; - partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi); + partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi); fraclo = comm->ysplit[comm->myloc[1]]; frachi = comm->ysplit[comm->myloc[1]+1]; - partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + partition_grid(ny,fraclo,frachi,shift_grid,0,inylo,inyhi); fraclo = comm->zsplit[comm->myloc[2]]; frachi = comm->zsplit[comm->myloc[2]+1]; - partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi); + partition_grid(nz,fraclo,frachi,shift_grid,zextra,inzlo,inzhi); } else { fraclo = comm->mysplit[0][0]; frachi = comm->mysplit[0][1]; - partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi); + partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi); fraclo = comm->mysplit[1][0]; frachi = comm->mysplit[1][1]; - partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi); + partition_grid(ny,fraclo,frachi,shift_grid,0,inylo,inyhi); fraclo = comm->mysplit[2][0]; frachi = comm->mysplit[2][1]; - partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi); + partition_grid(nz,fraclo,frachi,shift_grid,zextra,inzlo,inzhi); } // extend owned grid bounds with ghost grid cells in each direction ghost_grid(); + // additional intialization + // other constructor invokes this directly + + initialize(); + + // 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; +} + +/* ---------------------------------------------------------------------- + additional one-time setup common to both constructors + ---------------------------------------------------------------------- */ + +void Grid3d::initialize() +{ // error check on size of grid stored by this proc bigint total = (bigint) @@ -347,7 +427,7 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, // default = caller grid is allocated to ghost grid // used when computing pack/unpack lists in indices() - // these values can be overridden by caller using set_larger_grid() + // these values can be overridden using set_caller_grid() fullxlo = outxlo; fullxhi = outxhi; @@ -381,22 +461,6 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, // 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; } /* ---------------------------------------------------------------------- @@ -406,11 +470,10 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, 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 - // NOTE: for future support of zfactor - 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 + 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 @@ -420,12 +483,19 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, ------------------------------------------------------------------------- */ void Grid3d::partition_grid(int ngrid, double fraclo, double frachi, - double shift, int &lo, int &hi) + double shift, int extra, int &lo, int &hi) { - lo = static_cast (fraclo * ngrid); - while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; - hi = static_cast (frachi * ngrid); - while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + if (extra == 0) { + lo = static_cast (fraclo * ngrid); + while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + hi = static_cast (frachi * ngrid); + while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + } else { + lo = static_cast (fraclo * ngrid/zfactor); + while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++; + hi = static_cast (frachi * ngrid/zfactor); + while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--; + } } /* ---------------------------------------------------------------------- @@ -468,6 +538,8 @@ void Grid3d::ghost_grid() double dxinv = nx / prd[0]; double dyinv = ny / prd[1]; double dzinv = nz / prd[2]; + double dzinv_slab = nz / (prd[2] * zfactor); + int lo, hi; lo = static_cast((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom_lo + OFFSET) - OFFSET; @@ -486,10 +558,33 @@ void Grid3d::ghost_grid() outzlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo); outzhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi); } else { + lo = static_cast((sublo[2]-dist[2]-boxlo[2]) * dzinv_slab + shift_atom_lo + OFFSET) - OFFSET; + hi = static_cast((subhi[2]+dist[2]-boxlo[2]) * dzinv_slab + shift_atom_hi + OFFSET) - OFFSET; + outzlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo); + outzhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi); } - - // limit out xyz lo/hi indices to global grid for non-periodic dims + // if zextra = 1: + // adjust grid boundaries for processors at +z end, + // to include added empty grid cells between periodically repeating slabs + // in this case: + // want grid data forward communicated from +z proc to -z proc, but not vice versa + // want grid data reverse communicated from -z proc to +z proc, but not vice versa + // this is accomplished by inzhi = outzhi on +z end (no ghost cells) + // also insure no other procs use ghost cells beyond +z limit + + if (zextra) { + if (layout != Comm::LAYOUT_TILED) { + if (comm->myloc[2] == comm->procgrid[2]-1) inzhi = outzhi = nz - 1; + } else { + if (comm->mysplit[2][1] == 1.0) inzhi = outzhi = nz - 1; + } + outzhi = MIN(outzhi,nz-1); + } + + // limit out xyz lo/hi indices to global grid for non-periodic dims + // if zextra = 1 (e.g. PPPM), grid is still fully periodic + int *periodicity = domain->periodicity; if (!periodicity[0]) { @@ -506,45 +601,6 @@ void Grid3d::ghost_grid() outzlo = MAX(0,outzlo); outzhi = MIN(nz-1,outzhi); } - - - // NOTE: zperiod effects, also in owned grid cells ? - - /* - - // extent of zprd when 2d slab mode is selected - - double zprd_slab = zprd*slab_volfactor; - - // 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; - nhi = static_cast ((subhi[2]+dist[2]-boxlo[2]) * - nz_pppm/zprd_slab + shift) - OFFSET; - nzlo_out = nlo + nlower; - nzhi_out = nhi + nupper; - - // for slab PPPM, change the grid boundary for processors at +z end - // to include the empty volume between periodically repeating slabs - // for slab PPPM, want charge data communicated from -z proc to +z proc, - // but not vice versa, also want field data communicated from +z proc to - // -z proc, but not vice versa - // this is accomplished by nzhi_in = nzhi_out on +z end (no ghost cells) - // also insure no other procs use ghost cells beyond +z limit - // differnet logic for non-tiled vs tiled decomposition - - if (slabflag == 1) { - if (comm->layout != Comm::LAYOUT_TILED) { - if (comm->myloc[2] == comm->procgrid[2]-1) nzhi_in = nzhi_out = nz_pppm - 1; - } else { - if (comm->mysplit[2][1] == 1.0) nzhi_in = nzhi_out = nz_pppm - 1; - } - nzhi_out = MIN(nzhi_out,nz_pppm-1); - } - - */ - } /* ---------------------------------------------------------------------- @@ -1770,9 +1826,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_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]); + partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,0,obox[0],obox[1]); + partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,0,obox[2],obox[3]); + partition_grid(nz,zsplit[k],zsplit[k+1],shift_grid,zextra,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 39734a482f..8c7ed114cc 100644 --- a/src/grid3d.h +++ b/src/grid3d.h @@ -23,6 +23,8 @@ class Grid3d : protected Pointers { enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes 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(); void set_distance(double); @@ -31,7 +33,7 @@ class Grid3d : protected Pointers { void set_shift_grid(double); void set_shift_atom(double, double); void set_zfactor(double); - void set_larger_grid(int, int, int, int, int, int); + void set_caller_grid(int, int, int, int, int, int); void set_proc_neighs(int, int, int, int, int, int); int identical(Grid3d *); @@ -59,7 +61,7 @@ class Grid3d : protected Pointers { MPI_Comm gridcomm; // communicator for this class // usually world, but MSM calls with subset - // inputs from caller + // input from caller int nx, ny, nz; // size of global grid in all 3 dims double maxdist; // distance owned atoms can move outside subdomain @@ -230,13 +232,11 @@ class Grid3d : protected Pointers { // internal methods // ------------------------------------------- - void partition_grid(int, double, double, double, int &, int &); + void initialize(); + void partition_grid(int, double, double, double, int, int &, int &); void ghost_grid(); void extract_comm_info(); - void store(int, int, int, int, int, int, int, int, int, int, int, int, - int, int, int, int, int, int, int, int, int, int, int, int); - void setup_comm_brick(int &, int &); void setup_comm_tiled(int &, int &); int ghost_adjacent_brick();