ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

debugging of grid remap

Browse Source
This commit is contained in:
Steve Plimpton
2022-10-27 16:40:53 -06:00
parent 7bf4c8d54a
commit b6e29fd5d7
8 changed files with 275 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
src/EXTRA-FIX/fix_ttm_grid.cpp
View File

@ -43,8 +43,8 @@ static constexpr int CHUNK = 1024;
// OFFSET avoids outside-of-box atoms being rounded to grid pts incorrectly
// SHIFT = 0.0 assigns atoms to lower-left grid pt
// SHIFT = 0.5 assigns atoms to nearest grid pt
// use SHIFT = 0.0 for now since it allows fix ave/chunk
// to spatially average consistent with the TTM grid
// use SHIFT = 0.0 to match fix ttm
// also it allows fix ave/chunk to spatially average consistently
static constexpr int OFFSET = 16384;
static constexpr double SHIFT = 0.5;
@ -61,6 +61,7 @@ FixTTMGrid::FixTTMGrid(LAMMPS *lmp, int narg, char **arg) :
if (outfile) error->all(FLERR,"Fix ttm/grid does not support outfile option - "
"use dump grid command or restart files instead");
shift = OFFSET + SHIFT;
skin_original = neighbor->skin;
}
@ -113,7 +114,7 @@ void FixTTMGrid::init()
FixTTM::init();
if (neighbor->skin > skin_original)
error->all(FLERR,"Cannot extend neighbor skin after fix ttm/griddefined");
error->all(FLERR,"Cannot extend neighbor skin after fix ttm/grid defined");
}
/* ---------------------------------------------------------------------- */
@ -499,6 +500,9 @@ void FixTTMGrid::reset_grid()
return;
} else delete gridnew;
// DEBUG
if (comm->me == 0) printf("Remapping grid on step %ld\n",update->ntimestep);
// delete grid data which doesn't need to persist from previous to new decomp
memory->destroy(grid_buf1);
@ -510,6 +514,9 @@ void FixTTMGrid::reset_grid()
grid_previous = grid;
T_electron_previous = T_electron;
nxlo_out_previous = nxlo_out;
nylo_out_previous = nylo_out;
nzlo_out_previous = nzlo_out;
// allocate new per-grid data for new decomposition
@ -530,12 +537,23 @@ void FixTTMGrid::reset_grid()
memory->destroy(remap_buf2);
// delete grid data and grid for previous decomposition
// NOTE: need to set offsets
int nxlo_out_prev,nylo_out_prev,nzlo_out_prev;
memory->destroy3d_offset(T_electron_previous,
nzlo_out_prev, nylo_out_prev, nxlo_out_prev);
nzlo_out_previous, nylo_out_previous,
nxlo_out_previous);
delete grid_previous;
// communicate temperatures to ghost cells on new grid
grid->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
// zero new net_energy_transfer
// in case compute_vector accesses it on timestep 0
outflag = 0;
memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
}
/* ----------------------------------------------------------------------
@ -593,7 +611,8 @@ void FixTTMGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *l
void FixTTMGrid::pack_remap_grid(void *vbuf, int nlist, int *list)
{
auto buf = (double *) vbuf;
double *src = &T_electron_previous[nzlo_out][nylo_out][nxlo_out];
double *src =
&T_electron_previous[nzlo_out_previous][nylo_out_previous][nxlo_out_previous];
for (int i = 0; i < nlist; i++) buf[i] = src[list[i]];
}

2
src/EXTRA-FIX/fix_ttm_grid.h
View File

@ -63,6 +63,7 @@ class FixTTMGrid : public FixTTM {
int nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out;
double delxinv, delyinv, delzinv;
double skin_original;
double shift;
FILE *fpout;
class Grid3d *grid;
@ -72,6 +73,7 @@ class FixTTMGrid : public FixTTM {
double *grid_buf1, *grid_buf2;
double ***T_electron_read;
int nxlo_out_previous,nylo_out_previous,nzlo_out_previous;
void allocate_grid() override;
void deallocate_grid() override;

2
src/comm.cpp
View File

@ -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

2
src/compute.cpp
View File

@ -62,7 +62,7 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) :
// set child class defaults
scalar_flag = vector_flag = array_flag = 0;
peratom_flag = local_flag = 0;
peratom_flag = local_flag = pergrid_flag = 0;
size_vector_variable = size_array_rows_variable = 0;
tempflag = pressflag = peflag = 0;

153
src/grid2d.cpp
View File

@ -298,7 +298,9 @@ Grid2d::~Grid2d()
memory->destroy(copy[i].unpacklist);
}
memory->sfree(copy);
delete [] requests;
delete [] requests_remap;
memory->sfree(rcbinfo);
@ -345,6 +347,7 @@ void Grid2d::store(int ixlo, int ixhi, int iylo, int iyhi,
recv = nullptr;
copy = nullptr;
requests = nullptr;
requests_remap = nullptr;
xsplit = ysplit = zsplit = nullptr;
grid2proc = nullptr;
@ -695,8 +698,8 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2)
double xlo,xhi,ylo,yhi;
int ghostbox[4],pbc[2];
// find overlaps of my extended ghost box with all other procs
// accounts for crossings of periodic boundaries
// find overlaps of my extended ghost box with all owned boxes
// accounts for ghost box overlapping periodic boundaries
// noverlap = # of overlaps, including self
// overlap = vector of overlap info using Overlap data struct
@ -708,11 +711,12 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2)
pbc[0] = pbc[1] = 0;
Overlap *overlap;
int noverlap = compute_overlap(ghostbox,pbc,overlap);
int noverlap = compute_overlap(1,ghostbox,pbc,overlap);
// send each proc an overlap message
// content: me, index of my overlap, box that overlaps with its owned cells
// ncopy = # of overlaps with myself, across a periodic boundary
// skip copy to self when non-PBC
int *proclist;
memory->create(proclist,noverlap,"grid2d:proclist");
@ -723,8 +727,10 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2)
ncopy = 0;
for (m = 0; m < noverlap; m++) {
if (overlap[m].proc == me) ncopy++;
else {
if (overlap[m].proc == me) {
if (overlap[m].pbc[0] == 0 && overlap[m].pbc[1] == 0) continue;
ncopy++;
} else {
proclist[nsend_request] = overlap[m].proc;
srequest[nsend_request].sender = me;
srequest[nsend_request].index = m;
@ -799,12 +805,14 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2)
nrecv = nrecv_response;
// create Copy data struct from overlaps with self
// skip copy to self when non-PBC
copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"grid2d:copy");
ncopy = 0;
for (m = 0; m < noverlap; m++) {
if (overlap[m].proc != me) continue;
if (overlap[m].pbc[0] == 0 && overlap[m].pbc[1] == 0) continue;
xlo = overlap[m].box[0];
xhi = overlap[m].box[1];
ylo = overlap[m].box[2];
@ -952,7 +960,7 @@ void Grid2d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which,
}
/* ----------------------------------------------------------------------
forward comm on regular grid of procs via list of swaps with 6 neighbor procs
forward comm for brick decomp via list of swaps with 6 neighbor procs
------------------------------------------------------------------------- */
template < class T >
@ -982,7 +990,7 @@ forward_comm_brick(T *ptr, int nper, int /*nbyte*/, int which,
}
/* ----------------------------------------------------------------------
forward comm on tiled grid decomp via Send/Recv lists of each neighbor proc
forward comm for tiled decomp via Send/Recv lists of each neighbor proc
------------------------------------------------------------------------- */
template < class T >
@ -1057,7 +1065,7 @@ void Grid2d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which,
}
/* ----------------------------------------------------------------------
reverse comm on regular grid of procs via list of swaps with 6 neighbor procs
reverse comm for brick decomp via list of swaps with 6 neighbor procs
------------------------------------------------------------------------- */
template < class T >
@ -1087,7 +1095,7 @@ reverse_comm_brick(T *ptr, int nper, int /*nbyte*/, int which,
}
/* ----------------------------------------------------------------------
reverse comm on tiled grid decomp via Send/Recv lists of each neighbor proc
reverse comm for tiled decomp via Send/Recv lists of each neighbor proc
------------------------------------------------------------------------- */
template < class T >
@ -1147,8 +1155,6 @@ 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)
{
int m;
@ -1159,7 +1165,11 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
deallocate_remap();
// compute overlaps of old decomp owned box with all owned boxes in new decomp
// set layout to current Comm layout
layout = comm->layout;
// overlaps of my 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
@ -1168,7 +1178,7 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
pbc[0] = pbc[1] = 0;
Overlap *overlap_old;
int noverlap_old = compute_overlap(oldbox,pbc,overlap_old);
int noverlap_old = compute_overlap(0,oldbox,pbc,overlap_old);
// use overlap_old to construct send and copy lists
@ -1193,11 +1203,11 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
send_remap[nsend_remap].npack =
old->indices(send_remap[nsend_remap].packlist,
box[0],box[1],box[2],box[3]);
}
nsend_remap++;
}
}
// compute overlaps of new decomp owned box with all owned boxes in old decomp
// overlaps of my new decomp owned box with all owned boxes in old decomp
// noverlap_new = # of overlaps, including self
// overlap_new = vector of overlap info in Overlap data struct
@ -1206,7 +1216,7 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
pbc[0] = pbc[1] = 0;
Overlap *overlap_new;
int noverlap_new = old->compute_overlap(newbox,pbc,overlap_new);
int noverlap_new = old->compute_overlap(0,newbox,pbc,overlap_new);
// use overlap_new to construct recv and copy lists
// set offsets for Recv data
@ -1228,9 +1238,9 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
recv_remap[nrecv_remap].nunpack =
indices(recv_remap[nrecv_remap].unpacklist,
box[0],box[1],box[2],box[3]);
}
nrecv_remap++;
}
}
// set offsets for received data
@ -1240,6 +1250,11 @@ void Grid2d::setup_remap(Grid2d *old, int &nremap_buf1, int &nremap_buf2)
offset += recv_remap[m].nunpack;
}
// length of MPI requests vector = nrecv_remap
delete [] requests_remap;
requests_remap = new MPI_Request[nrecv_remap];
// clean-up
clean_overlap();
@ -1348,11 +1363,13 @@ void Grid2d::read_file_style(T *ptr, FILE *fp, int nchunk, int maxline)
while (nread < ntotal) {
int nchunk = MIN(ntotal - nread, nchunk);
int eof = utils::read_lines_from_file(fp, nchunk, maxline, buffer, comm->me, world);
int eof = utils::read_lines_from_file(fp, nchunk, maxline, buffer, me, world);
if (eof) error->all(FLERR, "Unexpected end of grid data file");
nread += ptr->unpack_read_grid(buffer);
}
delete [] buffer;
}
/* ----------------------------------------------------------------------
@ -1377,8 +1394,6 @@ template < class T >
void Grid2d::write_file_style(T *ptr, int which,
int nper, int nbyte, MPI_Datatype datatype)
{
int me = comm->me;
// maxsize = max size of grid data owned by any proc
int mysize = (inxhi-inxlo+1) * (inyhi-inylo+1);
@ -1445,6 +1460,9 @@ void Grid2d::write_file_style(T *ptr, int which,
/* ----------------------------------------------------------------------
compute list of overlaps between box and the owned grid boxes of all procs
ghostflag = 1 if box includes ghost grid pts, called by setup_tiled()
ghostflag = 0 if box has no ghost grid pts, called by setup_remap()
layout != LAYOUT_TILED is only invoked by setup_remap()
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)
@ -1452,12 +1470,19 @@ void Grid2d::write_file_style(T *ptr, int which,
pbc = flags for grid periodicity in each dim
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 # of overlaps including with self, caller handles self overlaps as needed
return list of overlaps
for setup_tiled() this is what box_drop() computes
entire box for each overlap
caller will determine extent of overlap using PBC info
for setup_remap(), return extent of overlap (no PBC info involved)
use proc_box_uniform() or tiled() and MAX/MIN to determine this
------------------------------------------------------------------------- */
int Grid2d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap)
{
int obox[4];
memory->create(overlap_procs,nprocs,"grid2d:overlap_procs");
noverlap_list = maxoverlap_list = 0;
overlap_list = nullptr;
@ -1466,16 +1491,16 @@ int Grid2d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
// 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 iproclo = proc_index_uniform(box[0],ngrid[0],0,comm->xsplit);
int iprochi = proc_index_uniform(box[1],ngrid[0],0,comm->xsplit);
int jproclo = proc_index_uniform(box[2],ngrid[1],1,comm->ysplit);
int jprochi = proc_index_uniform(box[3],ngrid[1],1,comm->ysplit);
int obox[4];
// compute extent of overlap of box with with each proc's obox
for (int j = jproclo; j <= jprochi; j++)
for (int i = iproclo; i <= iprochi; i++) {
find_proc_box(i,j,obox);
proc_box_uniform(i,j,obox);
if (noverlap_list == maxoverlap_list) grow_overlap();
overlap_list[noverlap_list].proc = grid2proc[i][j][0];
@ -1486,8 +1511,26 @@ int Grid2d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
noverlap_list++;
}
} else if (layout == Comm::LAYOUT_TILED) {
} else {
box_drop(box,pbc);
// compute extent of overlap of box with with each proc's obox
if (ghostflag == 0) {
for (int m = 0; m < noverlap_list; m++) {
obox[0] = 0;
obox[1] = ngrid[0]-1;
obox[2] = 0;
obox[3] = ngrid[1]-1;
proc_box_tiled(overlap_list[m].proc,0,nprocs-1,obox);
overlap_list[m].box[0] = MAX(box[0],obox[0]);
overlap_list[m].box[1] = MIN(box[1],obox[1]);
overlap_list[m].box[2] = MAX(box[2],obox[2]);
overlap_list[m].box[3] = MIN(box[3],obox[3]);
}
}
}
overlap = overlap_list;
@ -1553,10 +1596,9 @@ void Grid2d::box_drop(int *box, int *pbc)
newpbc[1]++;
// box is not split, drop on RCB tree
// returns nprocs = # of procs it overlaps, including self
// returns np = # of procs it overlaps, including self
// returns proc_overlap = list of proc IDs it overlaps
// skip self overlap if no crossing of periodic boundaries
// do not skip self if overlap is in another periodic image
// add each overlap to overlap list
} else {
splitflag = 0;
@ -1564,8 +1606,6 @@ void Grid2d::box_drop(int *box, int *pbc)
box_drop_grid(box,0,nprocs-1,np,overlap_procs);
for (m = 0; m < np; m++) {
if (noverlap_list == maxoverlap_list) grow_overlap();
if (overlap_procs[m] == me &&
pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0) continue;
overlap_list[noverlap_list].proc = overlap_procs[m];
for (i = 0; i < 4; i++) overlap_list[noverlap_list].box[i] = box[i];
for (i = 0; i < 2; i++) overlap_list[noverlap_list].pbc[i] = pbc[i];
@ -1648,11 +1688,11 @@ void Grid2d::deallocate_remap()
{
for (int i = 0; i < nsend_remap; i++)
memory->destroy(send_remap[i].packlist);
memory->sfree(send_remap);
delete [] send_remap;
for (int i = 0; i < nrecv_remap; i++)
memory->destroy(recv_remap[i].unpacklist);
memory->sfree(recv_remap);
delete [] recv_remap;
if (self_remap) {
memory->destroy(copy_remap.packlist);
@ -1684,14 +1724,14 @@ int Grid2d::indices(int *&list, int xlo, int xhi, int ylo, int yhi)
}
/* ----------------------------------------------------------------------
find the comm->procgrid index for which proc owns the igrid index
find the comm->procgrid index = which proc owns the igrid index
igrid = grid index (0 to N-1) in dim
n = # of grid points in dim
dim = which dimension (0,1)
split = comm->x/y/z split for fractional bounds of each proc domain
------------------------------------------------------------------------- */
int Grid2d::find_proc_index(int igrid, int n, int dim, double *split)
int Grid2d::proc_index_uniform(int igrid, int n, int dim, double *split)
{
int gridlo,gridhi;
double fraclo,frachi;
@ -1716,13 +1756,13 @@ int Grid2d::find_proc_index(int igrid, int n, int dim, double *split)
}
/* ----------------------------------------------------------------------
find the grid box for proc with grid indices i,j
compute 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)
void Grid2d::proc_box_uniform(int i, int j, int *box)
{
int lo,hi;
double fraclo,frachi;
@ -1745,3 +1785,36 @@ void Grid2d::find_proc_box(int i, int j, int *box)
box[2] = lo;
box[3] = hi;
}
/* ----------------------------------------------------------------------
compute the grid box for proc within tiled decomposition
performed recursively until proclower = procupper = proc
return box = lo/hi bounds of proc's box in 2 dims
------------------------------------------------------------------------- */
void Grid2d::proc_box_tiled(int proc, int proclower, int procupper, int *box)
{
// end recursion when partition is a single proc
if (proclower == procupper) return;
// split processor partition
// procmid = 1st processor in upper half of partition
// = location in tree that stores this cut
// cut = index of first grid cell in upper partition
// dim = 0,1 dimension of cut
int procmid = proclower + (procupper - proclower) / 2 + 1;
int dim = rcbinfo[procmid].dim;
int cut = rcbinfo[procmid].cut;
// adjust box to reflect which half of partition the proc is in
if (proc < procmid) {
box[2*dim+1] = cut-1;
proc_box_tiled(proc,proclower,procmid-1,box);
} else {
box[2*dim] = cut;
proc_box_tiled(proc,procmid,procupper,box);
}
}

7
src/grid2d.h
View File

@ -226,7 +226,7 @@ protected:
template <class T> void read_file_style(T *, FILE *, int, int);
template <class T> void write_file_style(T *, int, int, int, MPI_Datatype);
int compute_overlap(int *, int *, Overlap *&);
int compute_overlap(int, int *, int *, Overlap *&);
void clean_overlap();
void box_drop(int *, int *);
void box_drop_grid(int *, int, int, int &, int *);
@ -236,8 +236,9 @@ protected:
void deallocate_remap();
int indices(int *&, int, int, int, int);
int find_proc_index(int, int, int, double *);
void find_proc_box(int, int, int *);
int proc_index_uniform(int, int, int, double *);
void proc_box_uniform(int, int, int *);
void proc_box_tiled(int, int, int, int *);
};
} // namespace LAMMPS_NS

155
src/grid3d.cpp
View File

@ -26,8 +26,6 @@
using namespace LAMMPS_NS;
enum{UNIFORM,TILED};
#define DELTA 16
static constexpr int OFFSET = 16384;
@ -317,7 +315,9 @@ Grid3d::~Grid3d()
memory->destroy(copy[i].unpacklist);
}
memory->sfree(copy);
delete [] requests;
delete [] requests_remap;
memory->sfree(rcbinfo);
@ -374,6 +374,7 @@ void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi,
recv = nullptr;
copy = nullptr;
requests = nullptr;
requests_remap = nullptr;
xsplit = ysplit = zsplit = nullptr;
grid2proc = nullptr;
@ -841,11 +842,12 @@ void Grid3d::setup_tiled(int &nbuf1, int &nbuf2)
pbc[0] = pbc[1] = pbc[2] = 0;
Overlap *overlap;
int noverlap = compute_overlap(ghostbox,pbc,overlap);
int noverlap = compute_overlap(1,ghostbox,pbc,overlap);
// send each proc an overlap message
// content: me, index of my overlap, box that overlaps with its owned cells
// ncopy = # of overlaps with myself, across a periodic boundary
// ncopy = # of overlaps with myself across a periodic boundary
// skip copy to self when non-PBC
int *proclist;
memory->create(proclist,noverlap,"grid3d:proclist");
@ -856,8 +858,11 @@ void Grid3d::setup_tiled(int &nbuf1, int &nbuf2)
ncopy = 0;
for (m = 0; m < noverlap; m++) {
if (overlap[m].proc == me) ncopy++;
else {
if (overlap[m].proc == me) {
if (overlap[m].pbc[0] == 0 && overlap[m].pbc[1] == 0 &&
overlap[m].pbc[2] == 0) continue;
ncopy++;
} else {
proclist[nsend_request] = overlap[m].proc;
srequest[nsend_request].sender = me;
srequest[nsend_request].index = m;
@ -939,12 +944,15 @@ void Grid3d::setup_tiled(int &nbuf1, int &nbuf2)
nrecv = nrecv_response;
// create Copy data struct from overlaps with self
// skip copy to self when non-PBC
copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"grid3d:copy");
ncopy = 0;
for (m = 0; m < noverlap; m++) {
if (overlap[m].proc != me) continue;
if (overlap[m].pbc[0] == 0 && overlap[m].pbc[1] == 0 &&
overlap[m].pbc[2] == 0) continue;
xlo = overlap[m].box[0];
xhi = overlap[m].box[1];
ylo = overlap[m].box[2];
@ -1303,7 +1311,11 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
deallocate_remap();
// compute overlaps of old decomp owned box with all owned boxes in new decomp
// set layout to current Comm layout
layout = comm->layout;
// overlaps of my 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
@ -1312,7 +1324,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
pbc[0] = pbc[1] = pbc[2] = 0;
Overlap *overlap_old;
int noverlap_old = compute_overlap(oldbox,pbc,overlap_old);
int noverlap_old = compute_overlap(0,oldbox,pbc,overlap_old);
// use overlap_old to construct send and copy lists
@ -1320,7 +1332,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
nsend_remap = 0;
for (m = 0; m < noverlap_old; m++) {
if (overlap_old[m].proc == me) self_remap =1;
if (overlap_old[m].proc == me) self_remap = 1;
else nsend_remap++;
}
@ -1338,11 +1350,11 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
send_remap[nsend_remap].npack =
old->indices(send_remap[nsend_remap].packlist,
box[0],box[1],box[2],box[3],box[4],box[5]);
}
nsend_remap++;
}
}
// compute overlaps of new decomp owned box with all owned boxes in old decomp
// overlaps of my new decomp owned box with all owned boxes in old decomp
// noverlap_new = # of overlaps, including self
// overlap_new = vector of overlap info in Overlap data struct
@ -1351,7 +1363,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
pbc[0] = pbc[1] = pbc[2] = 0;
Overlap *overlap_new;
int noverlap_new = old->compute_overlap(newbox,pbc,overlap_new);
int noverlap_new = old->compute_overlap(0,newbox,pbc,overlap_new);
// use overlap_new to construct recv and copy lists
// set offsets for Recv data
@ -1374,18 +1386,23 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
recv_remap[nrecv_remap].nunpack =
indices(recv_remap[nrecv_remap].unpacklist,
box[0],box[1],box[2],box[3],box[4],box[5]);
}
nrecv_remap++;
}
}
// set offsets for received data
int offset = 0;
for (m = 0; m < nrecv_remap; m++) {
recv[m].offset = offset;
recv_remap[m].offset = offset;
offset += recv_remap[m].nunpack;
}
// length of MPI requests vector = nrecv_remap
delete [] requests_remap;
requests_remap = new MPI_Request[nrecv_remap];
// clean-up
clean_overlap();
@ -1494,7 +1511,7 @@ void Grid3d::read_file_style(T *ptr, FILE *fp, int nchunk, int maxline)
while (nread < ntotal) {
int nchunk = MIN(ntotal - nread, nchunk);
int eof = utils::read_lines_from_file(fp, nchunk, maxline, buffer, comm->me, world);
int eof = utils::read_lines_from_file(fp, nchunk, maxline, buffer, me, world);
if (eof) error->all(FLERR, "Unexpected end of grid data file");
nread += ptr->unpack_read_grid(buffer);
@ -1525,8 +1542,6 @@ template < class T >
void Grid3d::write_file_style(T *ptr, int which,
int nper, int nbyte, MPI_Datatype datatype)
{
int me = comm->me;
// maxsize = max size of grid data owned by any proc
int mysize = (inxhi-inxlo+1) * (inyhi-inylo+1) * (inzhi-inzlo+1);
@ -1597,6 +1612,9 @@ void Grid3d::write_file_style(T *ptr, int which,
/* ----------------------------------------------------------------------
compute list of overlaps between box and the owned grid boxes of all procs
ghostflag = 1 if box includes ghost grid pts, called by setup_tiled()
ghostflag = 0 if box has no ghost grid pts, called by setup_remap()
layout != LAYOUT_TILED is only invoked by setup_remap()
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)
@ -1604,12 +1622,19 @@ void Grid3d::write_file_style(T *ptr, int which,
pbc = flags for grid periodicity in each dim
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 # of overlaps including with self, caller handles self overlaps as needed
return list of overlaps
for setup_tiled() this is what box_drop() computes
entire box for each overlap
caller will determine extent of overlap using PBC info
for setup_remap(), return extent of overlap (no PBC info involved)
use proc_box_uniform() or tiled() and MAX/MIN to determine this
------------------------------------------------------------------------- */
int Grid3d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap)
{
int obox[6];
memory->create(overlap_procs,nprocs,"grid3d:overlap_procs");
noverlap_list = maxoverlap_list = 0;
overlap_list = nullptr;
@ -1618,19 +1643,19 @@ int Grid3d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
// find comm->procgrid indices in each dim for box bounds
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 iproclo = proc_index_uniform(box[0],ngrid[0],0,xsplit);
int iprochi = proc_index_uniform(box[1],ngrid[0],0,xsplit);
int jproclo = proc_index_uniform(box[2],ngrid[1],1,ysplit);
int jprochi = proc_index_uniform(box[3],ngrid[1],1,ysplit);
int kproclo = proc_index_uniform(box[4],ngrid[2],2,zsplit);
int kprochi = proc_index_uniform(box[5],ngrid[2],2,zsplit);
int obox[6];
// compute extent of overlap of box with with each proc's obox
for (int k = kproclo; k <= kprochi; k++)
for (int j = jproclo; j <= jprochi; j++)
for (int i = iproclo; i <= iprochi; i++) {
find_proc_box(i,j,k,obox);
proc_box_uniform(i,j,k,obox);
if (noverlap_list == maxoverlap_list) grow_overlap();
overlap_list[noverlap_list].proc = grid2proc[i][j][k];
@ -1643,8 +1668,30 @@ int Grid3d::compute_overlap(int *box, int *pbc, Overlap *&overlap)
noverlap_list++;
}
} else if (layout == TILED) {
} else {
box_drop(box,pbc);
// compute extent of overlap of box with with each proc's obox
if (ghostflag == 0) {
for (int m = 0; m < noverlap_list; m++) {
obox[0] = 0;
obox[1] = ngrid[0]-1;
obox[2] = 0;
obox[3] = ngrid[1]-1;
obox[4] = 0;
obox[5] = ngrid[2]-1;
proc_box_tiled(overlap_list[m].proc,0,nprocs-1,obox);
overlap_list[m].box[0] = MAX(box[0],obox[0]);
overlap_list[m].box[1] = MIN(box[1],obox[1]);
overlap_list[m].box[2] = MAX(box[2],obox[2]);
overlap_list[m].box[3] = MIN(box[3],obox[3]);
overlap_list[m].box[4] = MAX(box[4],obox[4]);
overlap_list[m].box[5] = MIN(box[5],obox[5]);
}
}
}
overlap = overlap_list;
@ -1720,10 +1767,9 @@ void Grid3d::box_drop(int *box, int *pbc)
newpbc[2]++;
// box is not split, drop on RCB tree
// returns nprocs = # of procs it overlaps, including self
// returns np = # of procs it overlaps, including self
// returns proc_overlap = list of proc IDs it overlaps
// skip self overlap if no crossing of periodic boundaries
// do not skip self if overlap is in another periodic image
// add each overlap to overlap list
} else {
splitflag = 0;
@ -1731,8 +1777,6 @@ void Grid3d::box_drop(int *box, int *pbc)
box_drop_grid(box,0,nprocs-1,np,overlap_procs);
for (m = 0; m < np; m++) {
if (noverlap_list == maxoverlap_list) grow_overlap();
if (overlap_procs[m] == me &&
pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0) continue;
overlap_list[noverlap_list].proc = overlap_procs[m];
for (i = 0; i < 6; i++) overlap_list[noverlap_list].box[i] = box[i];
for (i = 0; i < 3; i++) overlap_list[noverlap_list].pbc[i] = pbc[i];
@ -1815,11 +1859,11 @@ void Grid3d::deallocate_remap()
{
for (int i = 0; i < nsend_remap; i++)
memory->destroy(send_remap[i].packlist);
memory->sfree(send_remap);
delete [] send_remap;
for (int i = 0; i < nrecv_remap; i++)
memory->destroy(recv_remap[i].unpacklist);
memory->sfree(recv_remap);
delete [] recv_remap;
if (self_remap) {
memory->destroy(copy_remap.packlist);
@ -1854,14 +1898,14 @@ int Grid3d::indices(int *&list,
}
/* ----------------------------------------------------------------------
find the comm->procgrid index for which proc owns the igrid index
find the comm->procgrid index = which proc owns the igrid index
igrid = grid index (0 to N-1) in dim
n = # of grid points in dim
dim = which dimension (0,1,2)
split = comm->x/y/z split for fractional bounds of each proc domain
------------------------------------------------------------------------- */
int Grid3d::find_proc_index(int igrid, int n, int dim, double *split)
int Grid3d::proc_index_uniform(int igrid, int n, int dim, double *split)
{
int gridlo,gridhi;
double fraclo,frachi;
@ -1886,13 +1930,13 @@ int Grid3d::find_proc_index(int igrid, int n, int dim, double *split)
}
/* ----------------------------------------------------------------------
find the grid box for proc with grid indices i,j,k
compute 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)
void Grid3d::proc_box_uniform(int i, int j, int k, int *box)
{
int lo,hi;
double fraclo,frachi;
@ -1924,3 +1968,36 @@ void Grid3d::find_proc_box(int i, int j, int k, int *box)
box[4] = lo;
box[5] = hi;
}
/* ----------------------------------------------------------------------
compute the grid box for proc within tiled decomposition
performed recursively until proclower = procupper = proc
return box = lo/hi bounds of proc's box in 3 dims
------------------------------------------------------------------------- */
void Grid3d::proc_box_tiled(int proc, int proclower, int procupper, int *box)
{
// end recursion when partition is a single proc
if (proclower == procupper) return;
// split processor partition
// procmid = 1st processor in upper half of partition
// = location in tree that stores this cut
// cut = index of first grid cell in upper partition
// dim = 0,1,2 dimension of cut
int procmid = proclower + (procupper - proclower) / 2 + 1;
int dim = rcbinfo[procmid].dim;
int cut = rcbinfo[procmid].cut;
// adjust box to reflect which half of partition the proc is in
if (proc < procmid) {
box[2*dim+1] = cut-1;
proc_box_tiled(proc,proclower,procmid-1,box);
} else {
box[2*dim] = cut;
proc_box_tiled(proc,procmid,procupper,box);
}
}

7
src/grid3d.h
View File

@ -231,7 +231,7 @@ class Grid3d : protected Pointers {
template <class T> void read_file_style(T *, FILE *, int, int);
template <class T> void write_file_style(T *, int, int, int, MPI_Datatype);
int compute_overlap(int *, int *, Overlap *&);
int compute_overlap(int, int *, int *, Overlap *&);
void clean_overlap();
void box_drop(int *, int *);
void box_drop_grid(int *, int, int, int &, int *);
@ -241,8 +241,9 @@ class Grid3d : protected Pointers {
void deallocate_remap();
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 *);
int proc_index_uniform(int, int, int, double *);
void proc_box_uniform(int, int, int, int *);
void proc_box_tiled(int, int, int, int *);
};
} // namespace LAMMPS_NS