From 7bf4c8d54aa806b12e977ef7023f0be7b5c6e333 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 26 Oct 2022 15:14:06 -0600 Subject: [PATCH] more debugging --- src/EXTRA-FIX/fix_ttm.cpp | 9 +- src/EXTRA-FIX/fix_ttm_grid.cpp | 4 +- src/comm.cpp | 2 +- src/fix.cpp | 4 +- src/grid2d.cpp | 216 ++++++++++++++++++--------------- src/grid2d.h | 16 ++- src/grid3d.cpp | 208 ++++++++++++++++++------------- src/grid3d.h | 7 ++ src/input.cpp | 2 - 9 files changed, 264 insertions(+), 204 deletions(-) diff --git a/src/EXTRA-FIX/fix_ttm.cpp b/src/EXTRA-FIX/fix_ttm.cpp index d5a919a6a2..5ad68edfbe 100644 --- a/src/EXTRA-FIX/fix_ttm.cpp +++ b/src/EXTRA-FIX/fix_ttm.cpp @@ -226,13 +226,6 @@ void FixTTM::init() if (domain->triclinic) error->all(FLERR,"Cannot use fix ttm with triclinic box"); - // to allow this, would have to reset grid bounds dynamically - // for RCB balancing would have to reassign grid pts to procs - // and create a new Grid3d, and pass old GC data to new GC - - if (domain->box_change) - error->all(FLERR,"Cannot use fix ttm with changing box shape, size, or sub-domains"); - // set force prefactors for (int i = 1; i <= atom->ntypes; i++) { @@ -323,7 +316,7 @@ void FixTTM::post_force(int /*vflag*/) flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5); flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5); flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - + f[i][0] += flangevin[i][0]; f[i][1] += flangevin[i][1]; f[i][2] += flangevin[i][2]; diff --git a/src/EXTRA-FIX/fix_ttm_grid.cpp b/src/EXTRA-FIX/fix_ttm_grid.cpp index 289fd675e0..240ebb9b3d 100644 --- a/src/EXTRA-FIX/fix_ttm_grid.cpp +++ b/src/EXTRA-FIX/fix_ttm_grid.cpp @@ -173,7 +173,7 @@ void FixTTMGrid::post_force(int /*vflag*/) f[i][2] += flangevin[i][2]; } } - + if (flag) error->one(FLERR,"Out of range fix ttm/grid atoms"); } @@ -498,7 +498,7 @@ void FixTTMGrid::reset_grid() delete gridnew; return; } else delete gridnew; - + // delete grid data which doesn't need to persist from previous to new decomp memory->destroy(grid_buf1); diff --git a/src/comm.cpp b/src/comm.cpp index 7e6c677c45..82f4ad0a9e 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -1,4 +1,4 @@ -// clang-format off +\// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories diff --git a/src/fix.cpp b/src/fix.cpp index 83e0650483..9b27486f48 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -81,8 +81,8 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : stores_ids = 0; scalar_flag = vector_flag = array_flag = 0; - peratom_flag = local_flag = 0; - global_freq = local_freq = peratom_freq = -1; + peratom_flag = local_flag = pergrid_flag = 0; + global_freq = local_freq = peratom_freq = pergrid_freq = -1; size_vector_variable = size_array_rows_variable = 0; comm_forward = comm_reverse = comm_border = 0; diff --git a/src/grid2d.cpp b/src/grid2d.cpp index f2ce0de0e7..54584b8008 100644 --- a/src/grid2d.cpp +++ b/src/grid2d.cpp @@ -26,8 +26,6 @@ using namespace LAMMPS_NS; -enum{BRICK,TILED}; - #define DELTA 16 static constexpr int OFFSET = 16384; @@ -73,9 +71,7 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, ny = gny; ngrid[0] = nx; ngrid[1] = ny; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // partition global grid across procs // i xyz lo/hi = lower/upper bounds of global grid this proc owns @@ -146,7 +142,7 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, // store grid bounds and proc neighs - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi, oxlo,oxhi,oylo,oyhi, @@ -189,13 +185,11 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, ny = gny; ngrid[0] = nx; ngrid[1] = ny; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // store grid bounds and proc neighs - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi, oxlo,oxhi,oylo,oyhi, @@ -240,14 +234,12 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag, ny = gny; ngrid[0] = nx; ngrid[1] = ny; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // store grid bounds and proc neighs if (flag == 1) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { // this assumes gcomm = world int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi, @@ -263,7 +255,7 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag, } } else if (flag == 2) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { store(ixlo,ixhi,iylo,iyhi, oxlo,oxhi,oylo,oyhi, oxlo,oxhi,oylo,oyhi, @@ -286,6 +278,11 @@ Grid2d::~Grid2d() } memory->sfree(swap); + delete [] xsplit; + delete [] ysplit; + delete [] zsplit; + memory->destroy(grid2proc); + // tiled comm data structs for (int i = 0; i < nsend; i++) @@ -338,31 +335,6 @@ void Grid2d::store(int ixlo, int ixhi, int iylo, int iyhi, fullylo = fylo; fullyhi = fyhi; - // for BRICK layout, proc xy lohi = my 4 neighbor procs in this MPI_Comm - - if (layout == BRICK) { - procxlo = pxlo; - procxhi = pxhi; - procylo = pylo; - procyhi = pyhi; - } - - // for TILED layout, create RCB tree of cut info for grid decomp - // access CommTiled to get cut dimension - // cut = this proc's inlo in that dim - // dim is -1 for proc 0, but never accessed - - if (layout == TILED) { - rcbinfo = (RCBinfo *) - memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); - RCBinfo rcbone; - rcbone.dim = comm->rcbcutdim; - if (rcbone.dim <= 0) rcbone.cut = inxlo; - else if (rcbone.dim == 1) rcbone.cut = inylo; - MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR, - rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm); - } - // internal data initializations nswap = maxswap = 0; @@ -373,8 +345,53 @@ void Grid2d::store(int ixlo, int ixhi, int iylo, int iyhi, recv = nullptr; copy = nullptr; requests = nullptr; - + + xsplit = ysplit = zsplit = nullptr; + grid2proc = nullptr; rcbinfo = nullptr; + + nsend_remap = nrecv_remap = self_remap = 0; + send_remap = nullptr; + recv_remap = nullptr; + + // for non TILED layout: + // proc xyz lohi = my 6 neighbor procs in this MPI_Comm + // xyz split = copy of 1d vectors in Comm + // grid2proc = copy of 3d array in Comm + + if (layout != Comm::LAYOUT_TILED) { + procxlo = pxlo; + procxhi = pxhi; + procylo = pylo; + procyhi = pyhi; + + xsplit = new double[comm->procgrid[0]+1]; + ysplit = new double[comm->procgrid[1]+1]; + memcpy(xsplit,comm->xsplit,(comm->procgrid[0]+1)*sizeof(double)); + memcpy(ysplit,comm->ysplit,(comm->procgrid[1]+1)*sizeof(double)); + + memory->create(grid2proc,comm->procgrid[0],comm->procgrid[1],comm->procgrid[2], + "grid3d:grid2proc"); + memcpy(&grid2proc[0][0][0],&comm->grid2proc[0][0][0], + comm->procgrid[0]*comm->procgrid[1]*comm->procgrid[2]*sizeof(int)); + } + + // for TILED layout: + // create RCB tree of cut info for grid decomp + // access CommTiled to get cut dimension + // cut = this proc's inlo in that dim + // dim is -1 for proc 0, but never accessed + + if (layout == Comm::LAYOUT_TILED) { + rcbinfo = (RCBinfo *) + memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); + RCBinfo rcbone; + rcbone.dim = comm->rcbcutdim; + if (rcbone.dim <= 0) rcbone.cut = inxlo; + else if (rcbone.dim == 1) rcbone.cut = inylo; + MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR, + rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm); + } } /* ---------------------------------------------------------------------- */ @@ -446,7 +463,7 @@ void Grid2d::get_bounds_ghost(int &xlo, int &xhi, int &ylo, int &yhi) void Grid2d::setup(int &nbuf1, int &nbuf2) { - if (layout == BRICK) setup_brick(nbuf1,nbuf2); + if (layout != Comm::LAYOUT_TILED) setup_brick(nbuf1,nbuf2); else setup_tiled(nbuf1,nbuf2); } @@ -865,7 +882,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) int Grid2d::ghost_adjacent() { - if (layout == BRICK) return ghost_adjacent_brick(); + if (layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); return ghost_adjacent_tiled(); } @@ -911,7 +928,7 @@ int Grid2d::ghost_adjacent_tiled() void Grid2d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { if (caller == KSPACE) forward_comm_brick((KSpace *) ptr,nper,nbyte,which, buf1,buf2,datatype); @@ -1016,7 +1033,7 @@ forward_comm_tiled(T *ptr, int nper, int nbyte, int which, void Grid2d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { if (caller == KSPACE) reverse_comm_brick((KSpace *) ptr,nper,nbyte,which, buf1,buf2,datatype); @@ -1130,28 +1147,18 @@ reverse_comm_tiled(T *ptr, int nper, int nbyte, int which, caller converts them to message size for grid data it stores ------------------------------------------------------------------------- */ +/* ------------------------------------------------------------------------- */ + void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2) -{ - if (layout == BRICK) setup_remap_brick(old,nremap_buf1,nremap_buf2); - else setup_remap_tiled(old,nremap_buf1,nremap_buf2); -} - -/* ------------------------------------------------------------------------- */ - -void Grid2d::setup_remap_brick(Grid2d *old, int &nremap_buf1, int &nremap_buf2) -{ - nremap_buf1 = 0; - nremap_buf2 = 0; -} - -/* ------------------------------------------------------------------------- */ - -void Grid2d::setup_remap_tiled(Grid2d *old, int &nremap_buf1, int &nremap_buf2) { int m; int pbc[2]; int *box; + // deallocated existing remap data structs + + deallocate_remap(); + // compute overlaps of old decomp owned box with all owned boxes in new decomp // noverlap_old = # of overlaps, including self // overlap_old = vector of overlap info in Overlap data struct @@ -1206,7 +1213,7 @@ void Grid2d::setup_remap_tiled(Grid2d *old, int &nremap_buf1, int &nremap_buf2) nrecv_remap = 0; for (m = 0; m < noverlap_new; m++) - if (overlap_old[m].proc != me) nrecv_remap++; + if (overlap_new[m].proc != me) nrecv_remap++; recv_remap = new Recv[nrecv_remap]; @@ -1438,12 +1445,12 @@ void Grid2d::write_file_style(T *ptr, int which, /* ---------------------------------------------------------------------- compute list of overlaps between box and the owned grid boxes of all procs - for brick decomp, done using Comm::grid2proc data struct - for tiled decomp, done via recursive box drop on RCB tree + for brick decomp of Grid, done using xyz split + grid2proc copied from Comm + for tiled decomp of Grid, done via recursive box drop on RCB tree box = 4 integers = (xlo,xhi,ylo,yhi) box can be owned cells or owned + ghost cells pbc = flags for grid periodicity in each dim - if box includes ghost cells, it can overlap PBCs + if box includes ghost cells, it can overlap PBCs (only for setup_tiled) each lo/hi value may extend beyond 0 to N-1 into another periodic image return # of overlaps including with self return list of overlaps @@ -1455,7 +1462,7 @@ int Grid2d::compute_overlap(int *box, int *pbc, Overlap *&overlap) noverlap_list = maxoverlap_list = 0; overlap_list = nullptr; - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { // find comm->procgrid indices in each dim for box bounds @@ -1464,42 +1471,22 @@ int Grid2d::compute_overlap(int *box, int *pbc, Overlap *&overlap) int jproclo = find_proc_index(box[2],ngrid[1],1,comm->ysplit); int jprochi = find_proc_index(box[3],ngrid[1],1,comm->ysplit); - // save comm->myloc values so can overwrite them k,j,i triple loop - // b/c comm->partition_grid uses comm->myloc + int obox[4]; - int save_myloc[3]; - save_myloc[0] = comm->myloc[0]; - save_myloc[1] = comm->myloc[1]; - save_myloc[2] = comm->myloc[2]; + for (int j = jproclo; j <= jprochi; j++) + for (int i = iproclo; i <= iprochi; i++) { + find_proc_box(i,j,obox); - int obox[6]; - - for (int k = 0; k <= 0; k++) - for (int j = jproclo; j <= jprochi; j++) - for (int i = iproclo; i <= iprochi; i++) { - comm->myloc[0] = i; - comm->myloc[1] = j; - comm->myloc[2] = k; - - comm->partition_grid(ngrid[0],ngrid[1],1,0.0, - obox[0],obox[1],obox[2],obox[3],obox[4],obox[5]); - - if (noverlap_list == maxoverlap_list) grow_overlap(); - overlap[noverlap_list].proc = comm->grid2proc[i][j][k]; - overlap[noverlap_list].box[0] = MAX(box[0],obox[0]); - overlap[noverlap_list].box[1] = MIN(box[1],obox[1]); - overlap[noverlap_list].box[2] = MAX(box[2],obox[2]); - overlap[noverlap_list].box[3] = MIN(box[3],obox[3]); - noverlap_list++; - } + if (noverlap_list == maxoverlap_list) grow_overlap(); + overlap_list[noverlap_list].proc = grid2proc[i][j][0]; + overlap_list[noverlap_list].box[0] = MAX(box[0],obox[0]); + overlap_list[noverlap_list].box[1] = MIN(box[1],obox[1]); + overlap_list[noverlap_list].box[2] = MAX(box[2],obox[2]); + overlap_list[noverlap_list].box[3] = MIN(box[3],obox[3]); + noverlap_list++; + } - // restore comm->myloc values - - comm->myloc[0] = save_myloc[0]; - comm->myloc[1] = save_myloc[1]; - comm->myloc[2] = save_myloc[2]; - - } else if (layout == TILED) { + } else if (layout == Comm::LAYOUT_TILED) { box_drop(box,pbc); } @@ -1727,3 +1714,34 @@ int Grid2d::find_proc_index(int igrid, int n, int dim, double *split) return m; } + +/* ---------------------------------------------------------------------- + find the grid box for proc with grid indices i,j + i,j,k = grid index (0 to N-1) in each dim + return lo/hi bounds of box in 2 dims + computation is same as Comm::partition_grid() +------------------------------------------------------------------------- */ + +void Grid2d::find_proc_box(int i, int j, int *box) +{ + int lo,hi; + double fraclo,frachi; + + fraclo = xsplit[i]; + frachi = xsplit[i+1]; + lo = static_cast (fraclo * ngrid[0]); + if (1.0*lo != fraclo*ngrid[0]) lo++; + hi = static_cast (frachi * ngrid[0]); + if (1.0*hi == frachi*ngrid[0]) hi--; + box[0] = lo; + box[1] = hi; + + fraclo = ysplit[j]; + frachi = ysplit[j+1]; + lo = static_cast (fraclo * ngrid[1]); + if (1.0*lo != fraclo*ngrid[1]) lo++; + hi = static_cast (frachi * ngrid[1]); + if (1.0*hi == frachi*ngrid[1]) hi--; + box[2] = lo; + box[3] = hi; +} diff --git a/src/grid2d.h b/src/grid2d.h index ca450a7610..15a4e245d3 100644 --- a/src/grid2d.h +++ b/src/grid2d.h @@ -173,6 +173,13 @@ protected: int *overlap_procs; // length of Nprocs in communicator + // BRICK decomposition + + double *xsplit,*ysplit,*zsplit; + int ***grid2proc; + + // TILED decomposition + // RCB tree of cut info // each proc contributes one value, except proc 0 @@ -187,10 +194,10 @@ protected: // includes overlaps across periodic boundaries, can also be self struct Overlap { - int proc; // proc whose owned cells overlap my ghost cells - int box[4]; // box that overlaps otherproc's owned cells + int proc; // proc whose cells overlap my cells + int box[4]; // box of my cells which overlap proc's cells // this box is wholly contained within global grid - int pbc[2]; // PBC offsets to convert box to a portion of my ghost box + int pbc[2]; // PBC offsets to convert my box to a portion of my ghost box // my ghost box may extend beyond global grid }; @@ -214,8 +221,6 @@ protected: template void reverse_comm_brick(T *, int, int, int, void *, void *, MPI_Datatype); template void reverse_comm_tiled(T *, int, int, int, void *, void *, MPI_Datatype); - void setup_remap_brick(Grid2d *, int &, int &); - void setup_remap_tiled(Grid2d *, int &, int &); template void remap_style(T *, int, int, void *, void *, MPI_Datatype); template void read_file_style(T *, FILE *, int, int); @@ -232,6 +237,7 @@ protected: int indices(int *&, int, int, int, int); int find_proc_index(int, int, int, double *); + void find_proc_box(int, int, int *); }; } // namespace LAMMPS_NS diff --git a/src/grid3d.cpp b/src/grid3d.cpp index aca0843457..2645ffae4a 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; -enum{BRICK,TILED}; +enum{UNIFORM,TILED}; #define DELTA 16 @@ -74,9 +74,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, nz = gnz; ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // partition global grid across procs // i xyz lo/hi = lower/upper bounds of global grid this proc owns @@ -158,7 +156,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, // store grid bounds and proc neighs - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi,izlo,izhi, oxlo,oxhi,oylo,oyhi,ozlo,ozhi, @@ -203,13 +201,11 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, nz = gnz; ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // store grid bounds and proc neighs - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi,izlo,izhi, oxlo,oxhi,oylo,oyhi,ozlo,ozhi, @@ -256,14 +252,12 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int flag, nz = gnz; ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz; - - if (comm->layout == Comm::LAYOUT_TILED) layout = TILED; - else layout = BRICK; + layout = comm->layout; // store grid bounds and proc neighs if (flag == 1) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { // this assumes gcomm = world int (*procneigh)[2] = comm->procneigh; store(ixlo,ixhi,iylo,iyhi,izlo,izhi, @@ -280,7 +274,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int flag, } } else if (flag == 2) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { store(ixlo,ixhi,iylo,iyhi,izlo,izhi, oxlo,oxhi,oylo,oyhi,ozlo,ozhi, oxlo,oxhi,oylo,oyhi,ozlo,ozhi, @@ -303,6 +297,11 @@ Grid3d::~Grid3d() } memory->sfree(swap); + delete [] xsplit; + delete [] ysplit; + delete [] zsplit; + memory->destroy(grid2proc); + // tiled comm data structs for (int i = 0; i < nsend; i++) @@ -365,34 +364,6 @@ void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi, fullzlo = fzlo; fullzhi = fzhi; - // for BRICK layout, proc xyz lohi = my 6 neighbor procs in this MPI_Comm - - if (layout == BRICK) { - procxlo = pxlo; - procxhi = pxhi; - procylo = pylo; - procyhi = pyhi; - proczlo = pzlo; - proczhi = pzhi; - } - - // for TILED layout, create RCB tree of cut info for grid decomp - // access CommTiled to get cut dimension - // cut = this proc's inlo in that dim - // dim is -1 for proc 0, but never accessed - - if (layout == TILED) { - rcbinfo = (RCBinfo *) - memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); - RCBinfo rcbone; - rcbone.dim = comm->rcbcutdim; - if (rcbone.dim <= 0) rcbone.cut = inxlo; - else if (rcbone.dim == 1) rcbone.cut = inylo; - else if (rcbone.dim == 2) rcbone.cut = inzlo; - MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR, - rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm); - } - // internal data initializations nswap = maxswap = 0; @@ -404,11 +375,57 @@ void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi, copy = nullptr; requests = nullptr; + xsplit = ysplit = zsplit = nullptr; + grid2proc = nullptr; rcbinfo = nullptr; nsend_remap = nrecv_remap = self_remap = 0; send_remap = nullptr; recv_remap = nullptr; + + // for non TILED layout: + // proc xyz lohi = my 6 neighbor procs in this MPI_Comm + // xyz split = copy of 1d vectors in Comm + // grid2proc = copy of 3d array in Comm + + if (layout != Comm::LAYOUT_TILED) { + procxlo = pxlo; + procxhi = pxhi; + procylo = pylo; + procyhi = pyhi; + proczlo = pzlo; + proczhi = pzhi; + + xsplit = new double[comm->procgrid[0]+1]; + ysplit = new double[comm->procgrid[1]+1]; + zsplit = new double[comm->procgrid[2]+1]; + memcpy(xsplit,comm->xsplit,(comm->procgrid[0]+1)*sizeof(double)); + memcpy(ysplit,comm->ysplit,(comm->procgrid[1]+1)*sizeof(double)); + memcpy(zsplit,comm->zsplit,(comm->procgrid[2]+1)*sizeof(double)); + + memory->create(grid2proc,comm->procgrid[0],comm->procgrid[1],comm->procgrid[2], + "grid3d:grid2proc"); + memcpy(&grid2proc[0][0][0],&comm->grid2proc[0][0][0], + comm->procgrid[0]*comm->procgrid[1]*comm->procgrid[2]*sizeof(int)); + } + + // for TILED layout: + // create RCB tree of cut info for grid decomp + // access CommTiled to get cut dimension + // cut = this proc's inlo in that dim + // dim is -1 for proc 0, but never accessed + + if (layout == Comm::LAYOUT_TILED) { + rcbinfo = (RCBinfo *) + memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); + RCBinfo rcbone; + rcbone.dim = comm->rcbcutdim; + if (rcbone.dim <= 0) rcbone.cut = inxlo; + else if (rcbone.dim == 1) rcbone.cut = inylo; + else if (rcbone.dim == 2) rcbone.cut = inzlo; + MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR, + rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm); + } } /* ---------------------------------------------------------------------- */ @@ -489,7 +506,7 @@ void Grid3d::get_bounds_ghost(int &xlo, int &xhi, int &ylo, int &yhi, void Grid3d::setup(int &nbuf1, int &nbuf2) { - if (layout == BRICK) setup_brick(nbuf1,nbuf2); + if (layout != Comm::LAYOUT_TILED) setup_brick(nbuf1,nbuf2); else setup_tiled(nbuf1,nbuf2); } @@ -1009,7 +1026,7 @@ void Grid3d::setup_tiled(int &nbuf1, int &nbuf2) int Grid3d::ghost_adjacent() { - if (layout == BRICK) return ghost_adjacent_brick(); + if (layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); return ghost_adjacent_tiled(); } @@ -1057,7 +1074,7 @@ int Grid3d::ghost_adjacent_tiled() void Grid3d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { if (caller == KSPACE) forward_comm_brick((KSpace *) ptr,nper,nbyte,which, buf1,buf2,datatype); @@ -1162,7 +1179,7 @@ forward_comm_tiled(T *ptr, int nper, int nbyte, int which, void Grid3d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { if (caller == KSPACE) reverse_comm_brick((KSpace *) ptr,nper,nbyte,which, buf1,buf2,datatype); @@ -1341,7 +1358,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2) nrecv_remap = 0; for (m = 0; m < noverlap_new; m++) - if (overlap_old[m].proc != me) nrecv_remap++; + if (overlap_new[m].proc != me) nrecv_remap++; recv_remap = new Recv[nrecv_remap]; @@ -1580,12 +1597,12 @@ void Grid3d::write_file_style(T *ptr, int which, /* ---------------------------------------------------------------------- compute list of overlaps between box and the owned grid boxes of all procs - for brick decomp, done using Comm::grid2proc data struct - for tiled decomp, done via recursive box drop on RCB tree + for brick decomp of Grid, done using xyz split + grid2proc copied from Comm + for tiled decomp of Grid, done via recursive box drop on RCB tree box = 6 integers = (xlo,xhi,ylo,yhi,zlo,zhi) box can be owned cells or owned + ghost cells pbc = flags for grid periodicity in each dim - if box includes ghost cells, it can overlap PBCs + if box includes ghost cells, it can overlap PBCs (only for setup_tiled) each lo/hi value may extend beyond 0 to N-1 into another periodic image return # of overlaps including with self return list of overlaps @@ -1597,54 +1614,35 @@ int Grid3d::compute_overlap(int *box, int *pbc, Overlap *&overlap) noverlap_list = maxoverlap_list = 0; overlap_list = nullptr; - if (layout == BRICK) { + if (layout != Comm::LAYOUT_TILED) { // find comm->procgrid indices in each dim for box bounds - int iproclo = find_proc_index(box[0],ngrid[0],0,comm->xsplit); - int iprochi = find_proc_index(box[1],ngrid[0],0,comm->xsplit); - int jproclo = find_proc_index(box[2],ngrid[1],1,comm->ysplit); - int jprochi = find_proc_index(box[3],ngrid[1],1,comm->ysplit); - int kproclo = find_proc_index(box[4],ngrid[2],2,comm->zsplit); - int kprochi = find_proc_index(box[5],ngrid[2],2,comm->zsplit); - - // save comm->myloc values so can overwrite them k,j,i triple loop - // b/c comm->partition_grid uses comm->myloc - - int save_myloc[3]; - save_myloc[0] = comm->myloc[0]; - save_myloc[1] = comm->myloc[1]; - save_myloc[2] = comm->myloc[2]; + int iproclo = find_proc_index(box[0],ngrid[0],0,xsplit); + int iprochi = find_proc_index(box[1],ngrid[0],0,xsplit); + int jproclo = find_proc_index(box[2],ngrid[1],1,ysplit); + int jprochi = find_proc_index(box[3],ngrid[1],1,ysplit); + int kproclo = find_proc_index(box[4],ngrid[2],2,zsplit); + int kprochi = find_proc_index(box[5],ngrid[2],2,zsplit); int obox[6]; for (int k = kproclo; k <= kprochi; k++) for (int j = jproclo; j <= jprochi; j++) for (int i = iproclo; i <= iprochi; i++) { - comm->myloc[0] = i; - comm->myloc[1] = j; - comm->myloc[2] = k; - - comm->partition_grid(ngrid[0],ngrid[1],ngrid[2],0.0, - obox[0],obox[1],obox[2],obox[3],obox[4],obox[5]); + find_proc_box(i,j,k,obox); if (noverlap_list == maxoverlap_list) grow_overlap(); - overlap[noverlap_list].proc = comm->grid2proc[i][j][k]; - overlap[noverlap_list].box[0] = MAX(box[0],obox[0]); - overlap[noverlap_list].box[1] = MIN(box[1],obox[1]); - overlap[noverlap_list].box[2] = MAX(box[2],obox[2]); - overlap[noverlap_list].box[3] = MIN(box[3],obox[3]); - overlap[noverlap_list].box[4] = MAX(box[4],obox[4]); - overlap[noverlap_list].box[5] = MIN(box[5],obox[5]); + overlap_list[noverlap_list].proc = grid2proc[i][j][k]; + overlap_list[noverlap_list].box[0] = MAX(box[0],obox[0]); + overlap_list[noverlap_list].box[1] = MIN(box[1],obox[1]); + overlap_list[noverlap_list].box[2] = MAX(box[2],obox[2]); + overlap_list[noverlap_list].box[3] = MIN(box[3],obox[3]); + overlap_list[noverlap_list].box[4] = MAX(box[4],obox[4]); + overlap_list[noverlap_list].box[5] = MIN(box[5],obox[5]); noverlap_list++; } - // restore comm->myloc values - - comm->myloc[0] = save_myloc[0]; - comm->myloc[1] = save_myloc[1]; - comm->myloc[2] = save_myloc[2]; - } else if (layout == TILED) { box_drop(box,pbc); } @@ -1759,7 +1757,7 @@ void Grid3d::box_drop(int *box, int *pbc) ------------------------------------------------------------------------- */ void Grid3d::box_drop_grid(int *box, int proclower, int procupper, - int &np, int *plist) + int &np, int *plist) { // end recursion when partition is a single proc // add proclower to plist @@ -1886,3 +1884,43 @@ int Grid3d::find_proc_index(int igrid, int n, int dim, double *split) return m; } + +/* ---------------------------------------------------------------------- + find the grid box for proc with grid indices i,j,k + i,j,k = grid index (0 to N-1) in each dim + return lo/hi bounds of box in 3 dims + computation is same as Comm::partition_grid() +------------------------------------------------------------------------- */ + +void Grid3d::find_proc_box(int i, int j, int k, int *box) +{ + int lo,hi; + double fraclo,frachi; + + fraclo = xsplit[i]; + frachi = xsplit[i+1]; + lo = static_cast (fraclo * ngrid[0]); + if (1.0*lo != fraclo*ngrid[0]) lo++; + hi = static_cast (frachi * ngrid[0]); + if (1.0*hi == frachi*ngrid[0]) hi--; + box[0] = lo; + box[1] = hi; + + fraclo = ysplit[j]; + frachi = ysplit[j+1]; + lo = static_cast (fraclo * ngrid[1]); + if (1.0*lo != fraclo*ngrid[1]) lo++; + hi = static_cast (frachi * ngrid[1]); + if (1.0*hi == frachi*ngrid[1]) hi--; + box[2] = lo; + box[3] = hi; + + fraclo = zsplit[k]; + frachi = zsplit[k+1]; + lo = static_cast (fraclo * ngrid[2]); + if (1.0*lo != fraclo*ngrid[2]) lo++; + hi = static_cast (frachi * ngrid[2]); + if (1.0*hi == frachi*ngrid[2]) hi--; + box[4] = lo; + box[5] = hi; +} diff --git a/src/grid3d.h b/src/grid3d.h index 3d653cd3b7..a9ed5be97f 100644 --- a/src/grid3d.h +++ b/src/grid3d.h @@ -179,6 +179,12 @@ class Grid3d : protected Pointers { int *overlap_procs; // length of Nprocs in communicator + // BRICK decomposition + + double *xsplit,*ysplit,*zsplit; + int ***grid2proc; + + // TILED decomposition // RCB tree of cut info // each proc contributes one value, except proc 0 @@ -236,6 +242,7 @@ class Grid3d : protected Pointers { int indices(int *&, int, int, int, int, int, int); int find_proc_index(int, int, int, double *); + void find_proc_box(int, int, int, int *); }; } // namespace LAMMPS_NS diff --git a/src/input.cpp b/src/input.cpp index 4c88579dfc..77571795fe 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -1438,10 +1438,8 @@ void Input::comm_style() } else if (strcmp(arg[0],"tiled") == 0) { if (comm->style == 1) return; Comm *oldcomm = comm; - if (lmp->kokkos) comm = new CommTiledKokkos(lmp,oldcomm); else comm = new CommTiled(lmp,oldcomm); - delete oldcomm; } else error->all(FLERR,"Illegal comm_style command"); }