From a541873b410733ca8e358671a6b246a87aabbfa9 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 21 Aug 2024 10:28:24 -0600 Subject: [PATCH] 2 bug fixes --- examples/grid/in.sph | 1 - src/fix_ave_grid.cpp | 59 +++++++++++++++++++++++++++++++------------- src/grid2d.cpp | 42 ++++++++++++++++--------------- src/grid2d.h | 3 +++ src/grid3d.cpp | 28 ++++++++++++++------- src/grid3d.h | 3 +++ 6 files changed, 89 insertions(+), 47 deletions(-) diff --git a/examples/grid/in.sph b/examples/grid/in.sph index 172b5be7d6..373aa50139 100644 --- a/examples/grid/in.sph +++ b/examples/grid/in.sph @@ -109,4 +109,3 @@ variable skin equal 0.3*${h} neighbor ${skin} bin run 6000 - diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index 16beed1123..53092900b3 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -684,11 +684,13 @@ void FixAveGrid::atom2grid() double ***vec3d = grid_sample->vec3d; double ****array3d = grid_sample->array3d; + // scan my owned atoms before tallying to bins // bin[i][dim] = indices of bin each atom is in + // count[][] or count[][][] = count of atoms contributing to each bin + // error check if any atom is out of bounds for my local+ghost grid cells // skip atom if group mask does not match - // check if any atom is out of bounds for my local grid - // for nonperiodic dim, remap atom to first/last bin if out of bounds - // count atoms contributing to each bin + // skip atom if out of bounds for a nonperiodic dim and discardflag = DISCARD + // if out of bounds for a nonperiodic dim and discardflag = KEEP, remap atom to first/last bin double *boxlo,*prd; int *periodicity = domain->periodicity; @@ -719,6 +721,7 @@ void FixAveGrid::atom2grid() if (triclinic) domain->x2lamda(nlocal); + // set to 1 int flag = 0; if (dimension == 2) { @@ -782,32 +785,54 @@ void FixAveGrid::atom2grid() iz = static_cast ((x[i][2]-boxlo[2])*dzinv + OFFSET) - OFFSET; if (ix < nxlo_out || ix > nxhi_out) { - if (periodicity[0]) flag = 1; - else if (discardflag == KEEP) { + if (periodicity[0]) { + flag = 1; + continue; + } else if (discardflag == KEEP) { if (ix < nxlo_out && nxlo_out == 0) ix = 0; else if (ix > nxhi_out && nxhi_out == nxgrid-1) ix = nxgrid-1; - else flag = 1; - } else skip[i] = 1; + else { + flag = 1; + continue; + } + } else { + skip[i] = 1; + continue; + } } if (iy < nylo_out || iy > nyhi_out) { - if (periodicity[1]) flag = 1; - else if (discardflag == KEEP) { + if (periodicity[1]) { + flag = 1; + continue; + } else if (discardflag == KEEP) { if (iy < nylo_out && nylo_out == 0) iy = 0; else if (iy > nyhi_out && nyhi_out == nygrid-1) iy = nygrid-1; - else flag = 1; - } else skip[i] = 1; - continue; + else { + flag = 1; + continue; + } + } else { + skip[i] = 1; + continue; + } } if (iz < nzlo_out || iz > nzhi_out) { - if (periodicity[2]) flag = 1; - else if (discardflag == KEEP) { + if (periodicity[2]) { + flag = 1; + continue; + } else if (discardflag == KEEP) { if (iz < nzlo_out && nzlo_out == 0) iz = 0; else if (iz > nzhi_out && nzhi_out == nzgrid-1) iz = nzgrid-1; - else flag = 1; - } else skip[i] = 1; - continue; + else { + flag = 1; + continue; + } + } else { + skip[i] = 1; + continue; + } } skip[i] = 0; diff --git a/src/grid2d.cpp b/src/grid2d.cpp index 4bac76d162..3b10bf036e 100644 --- a/src/grid2d.cpp +++ b/src/grid2d.cpp @@ -71,6 +71,11 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) : shift_atom_lo = shift_atom_hi = 0.0; yextra = 0; yfactor = 1.0; + + // layout_grid = how this grid instance is distributed across procs + // depends on comm->layout at time this Grid2d instance is created + + layout_grid = comm->layout; } /* ---------------------------------------------------------------------- @@ -112,6 +117,11 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int ixlo, int ixhi outylo = oylo; outyhi = oyhi; + // layout_grid = how this grid instance is distributed across procs + // depends on comm->layout at time this Grid2d instance is created + + layout_grid = comm->layout; + // additional intialization // other constructor invokes this from setup_grid() @@ -338,7 +348,7 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, int &oxlo, i double fraclo, frachi; - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { fraclo = comm->xsplit[comm->myloc[0]]; frachi = comm->xsplit[comm->myloc[0] + 1]; partition_grid(nx, fraclo, frachi, shift_grid, 0, inxlo, inxhi); @@ -417,8 +427,8 @@ void Grid2d::initialize() 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 + // two Grid instances will exist for duration of remap + // new Grid instance will then replace old one extract_comm_info(); } @@ -523,7 +533,7 @@ void Grid2d::ghost_grid() // also ensure no other procs use ghost cells beyond +y limit if (yextra) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != 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; @@ -560,7 +570,7 @@ void Grid2d::extract_comm_info() // xyz split = copy of 1d vectors in Comm // grid2proc = copy of 3d array in Comm - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { procxlo = comm->procneigh[0][0]; procxhi = comm->procneigh[0][1]; procylo = comm->procneigh[1][0]; @@ -584,7 +594,7 @@ void Grid2d::extract_comm_info() // RCBinfo.cut = this proc's inlo in that dim // Allgather creates the tree of dims and cuts - if (comm->layout == Comm::LAYOUT_TILED) { + if (layout_grid == Comm::LAYOUT_TILED) { rcbinfo = (RCBinfo *) memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); RCBinfo rcbone; @@ -593,9 +603,6 @@ void Grid2d::extract_comm_info() else if (rcbone.dim == 1) rcbone.cut = inylo; MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR, rcbinfo,sizeof(RCBinfo),MPI_CHAR,gridcomm); - printf("G2D ExtractCommInfo: 0 dim/cut: %d %d, 1 dim/cut: %d %d\n", - rcbinfo[0].dim,rcbinfo[0].cut, - rcbinfo[1].dim,rcbinfo[1].cut); } } @@ -617,7 +624,7 @@ void Grid2d::extract_comm_info() void Grid2d::setup_comm(int &nbuf1, int &nbuf2) { - if (comm->layout != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2); + if (layout_grid != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2); else setup_comm_tiled(nbuf1,nbuf2); } @@ -1041,7 +1048,7 @@ void Grid2d::setup_comm_tiled(int &nbuf1, int &nbuf2) int Grid2d::ghost_adjacent() { - if (comm->layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); + if (layout_grid != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); return ghost_adjacent_tiled(); } @@ -1087,7 +1094,7 @@ int Grid2d::ghost_adjacent_tiled() void Grid2d::forward_comm(int caller, void *ptr, int which, int nper, int nbyte, void *buf1, void *buf2, MPI_Datatype datatype) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { if (caller == KSPACE) forward_comm_brick((KSpace *) ptr,which,nper,nbyte, buf1,buf2,datatype); @@ -1192,7 +1199,7 @@ forward_comm_tiled(T *ptr, int which, int nper, int nbyte, void Grid2d::reverse_comm(int caller, void *ptr, int which, int nper, int nbyte, void *buf1, void *buf2, MPI_Datatype datatype) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { if (caller == KSPACE) reverse_comm_brick((KSpace *) ptr,which,nper,nbyte, buf1,buf2,datatype); @@ -1652,7 +1659,7 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap // test obox against appropriate layout - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { // find comm->procgrid indices in each dim for box bounds @@ -1803,7 +1810,6 @@ void Grid2d::box_drop_grid(int *box, int proclower, int procupper, if (proclower == procupper) { plist[np++] = proclower; - printf(" proc %d plist np %d plist %d\n",comm->me,np,plist[np-1]); return; } @@ -1815,13 +1821,9 @@ void Grid2d::box_drop_grid(int *box, int proclower, int procupper, // dim = 0,1,2 dimension of cut int procmid = proclower + (procupper - proclower) / 2 + 1; - - printf("G2D BoxDropGrid: proc %d proc lum: %d %d %d np %d\n",comm->me,proclower,procupper,procmid,np); - //printf(" dim %d cut %d\n",rcbinfo[procmid].dim,rcbinfo[procmid].cut); - int dim = rcbinfo[procmid].dim; int cut = rcbinfo[procmid].cut; - + if (box[2*dim] < cut) box_drop_grid(box,proclower,procmid-1,np,plist); if (box[2*dim+1] >= cut) box_drop_grid(box,procmid,procupper,np,plist); } diff --git a/src/grid2d.h b/src/grid2d.h index 8316f840be..9b8b3bb379 100644 --- a/src/grid2d.h +++ b/src/grid2d.h @@ -57,6 +57,9 @@ class Grid2d : protected Pointers { int me, nprocs; MPI_Comm gridcomm; // communicator for this class // usually world, but MSM calls with subset + int layout_grid; // how this grid instance is distributed across procs + // uses enum options for comm->layout + // load balancing can create a new Grid with new layout_grid // inputs from caller via constructor diff --git a/src/grid3d.cpp b/src/grid3d.cpp index f11e5d0513..e3f50a7b62 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -74,6 +74,11 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) : shift_atom_lo = shift_atom_hi = 0.0; zextra = 0; zfactor = 1.0; + + // layout_grid = how this grid instance is distributed across procs + // depends on comm->layout at time this Grid3d instance is created + + layout_grid = comm->layout; } /* ---------------------------------------------------------------------- @@ -122,6 +127,11 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz, outzlo = ozlo; outzhi = ozhi; + // layout_grid = how this grid instance is distributed across procs + // depends on comm->layout at time this Grid3d instance is created + + layout_grid = comm->layout; + // additional intialization // other constructor invokes this from setup_grid() @@ -375,7 +385,7 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi, double fraclo,frachi; - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { fraclo = comm->xsplit[comm->myloc[0]]; frachi = comm->xsplit[comm->myloc[0]+1]; partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi); @@ -579,7 +589,7 @@ void Grid3d::ghost_grid() // also ensure no other procs use ghost cells beyond +z limit if (zextra) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != 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; @@ -621,7 +631,7 @@ void Grid3d::extract_comm_info() // xyz split = copy of 1d vectors in Comm // grid2proc = copy of 3d array in Comm - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { procxlo = comm->procneigh[0][0]; procxhi = comm->procneigh[0][1]; procylo = comm->procneigh[1][0]; @@ -649,7 +659,7 @@ void Grid3d::extract_comm_info() // RCBinfo.cut = this proc's inlo in that dim // Allgather creates the tree of dims and cuts - if (comm->layout == Comm::LAYOUT_TILED) { + if (layout_grid == Comm::LAYOUT_TILED) { rcbinfo = (RCBinfo *) memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); RCBinfo rcbone; @@ -680,7 +690,7 @@ void Grid3d::extract_comm_info() void Grid3d::setup_comm(int &nbuf1, int &nbuf2) { - if (comm->layout != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2); + if (layout_grid != Comm::LAYOUT_TILED) setup_comm_brick(nbuf1,nbuf2); else setup_comm_tiled(nbuf1,nbuf2); } @@ -1207,7 +1217,7 @@ void Grid3d::setup_comm_tiled(int &nbuf1, int &nbuf2) int Grid3d::ghost_adjacent() { - if (comm->layout != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); + if (layout_grid != Comm::LAYOUT_TILED) return ghost_adjacent_brick(); return ghost_adjacent_tiled(); } @@ -1255,7 +1265,7 @@ int Grid3d::ghost_adjacent_tiled() void Grid3d::forward_comm(int caller, void *ptr, int which, int nper, int nbyte, void *buf1, void *buf2, MPI_Datatype datatype) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { if (caller == KSPACE) forward_comm_brick((KSpace *) ptr,which,nper,nbyte, buf1,buf2,datatype); @@ -1360,7 +1370,7 @@ forward_comm_tiled(T *ptr, int which, int nper, int nbyte, void Grid3d::reverse_comm(int caller, void *ptr, int which, int nper, int nbyte, void *buf1, void *buf2, MPI_Datatype datatype) { - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { if (caller == KSPACE) reverse_comm_brick((KSpace *) ptr,which,nper,nbyte, buf1,buf2,datatype); @@ -1825,7 +1835,7 @@ int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap return noverlap_list; } - if (comm->layout != Comm::LAYOUT_TILED) { + if (layout_grid != Comm::LAYOUT_TILED) { // find comm->procgrid indices in each dim for box bounds diff --git a/src/grid3d.h b/src/grid3d.h index 6a15c2c942..8a15b83d1f 100644 --- a/src/grid3d.h +++ b/src/grid3d.h @@ -59,6 +59,9 @@ class Grid3d : protected Pointers { int me, nprocs; MPI_Comm gridcomm; // communicator for this class // usually world, but MSM calls with subset + int layout_grid; // how this grid instance is distributed across procs + // uses enum options for comm->layout + // load balancing can create a new Grid with new layout_grid // input from caller