git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12752 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp
2014-11-20 18:02:47 +00:00
parent 973b24564e
commit 7f8f2fd3ed
9 changed files with 181 additions and 136 deletions

View File

@ -339,7 +339,7 @@ void Balance::command(int narg, char **arg)
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (domain->triclinic) domain->x2lamda(atom->nlocal);
Irregular *irregular = new Irregular(lmp); Irregular *irregular = new Irregular(lmp);
if (style == BISECTION) irregular->migrate_atoms(1,rcb->sendproc); if (style == BISECTION) irregular->migrate_atoms(1,1,rcb->sendproc);
else irregular->migrate_atoms(1); else irregular->migrate_atoms(1);
delete irregular; delete irregular;
if (domain->triclinic) domain->lamda2x(atom->nlocal); if (domain->triclinic) domain->lamda2x(atom->nlocal);
@ -537,7 +537,8 @@ int *Balance::bisection(int sortflag)
comm->rcbnew = 1; comm->rcbnew = 1;
int idim = rcb->cutdim; int idim = rcb->cutdim;
comm->rcbcutfrac = (rcb->cut - boxlo[idim]) / prd[idim]; if (idim >= 0) comm->rcbcutfrac = (rcb->cut - boxlo[idim]) / prd[idim];
else comm->rcbcutfrac = 0.0;
comm->rcbcutdim = idim; comm->rcbcutdim = idim;
double (*mysplit)[2] = comm->mysplit; double (*mysplit)[2] = comm->mysplit;

View File

@ -44,6 +44,7 @@ enum{SINGLE,MULTI}; // same as in Comm sub-styles
enum{MULTIPLE}; // same as in ProcMap enum{MULTIPLE}; // same as in ProcMap
enum{ONELEVEL,TWOLEVEL,NUMA,CUSTOM}; enum{ONELEVEL,TWOLEVEL,NUMA,CUSTOM};
enum{CART,CARTREORDER,XYZ}; enum{CART,CARTREORDER,XYZ};
enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -565,6 +566,84 @@ void Comm::set_proc_grid(int outflag)
} }
} }
/* ----------------------------------------------------------------------
determine which proc owns atom with coord x[3] based on current decomp
x will be in box (orthogonal) or lamda coords (triclinic)
if layout = UNIFORM, calculate owning proc directly
if layout = NONUNIFORM, iteratively find owning proc via binary search
if layout = TILED, CommTiled has its own method
return owning proc ID via grid2proc
return igx,igy,igz = logical grid loc of owing proc within 3d grid of procs
------------------------------------------------------------------------- */
int Comm::coord2proc(double *x, int &igx, int &igy, int &igz)
{
double *prd = domain->prd;
double *boxlo = domain->boxlo;
if (layout == LAYOUT_UNIFORM) {
if (triclinic == 0) {
igx = static_cast<int> (procgrid[0] * (x[0]-boxlo[0]) / prd[0]);
igy = static_cast<int> (procgrid[1] * (x[1]-boxlo[1]) / prd[1]);
igz = static_cast<int> (procgrid[2] * (x[2]-boxlo[2]) / prd[2]);
} else {
igx = static_cast<int> (procgrid[0] * x[0]);
igy = static_cast<int> (procgrid[1] * x[1]);
igz = static_cast<int> (procgrid[2] * x[2]);
}
} else if (layout == LAYOUT_NONUNIFORM) {
if (triclinic == 0) {
igx = binary((x[0]-boxlo[0])/prd[0],procgrid[0],xsplit);
igy = binary((x[1]-boxlo[1])/prd[1],procgrid[1],ysplit);
igz = binary((x[2]-boxlo[2])/prd[2],procgrid[2],zsplit);
} else {
igx = binary(x[0],procgrid[0],xsplit);
igy = binary(x[1],procgrid[1],ysplit);
igz = binary(x[2],procgrid[2],zsplit);
}
}
if (igx < 0) igx = 0;
if (igx >= procgrid[0]) igx = procgrid[0] - 1;
if (igy < 0) igy = 0;
if (igy >= procgrid[1]) igy = procgrid[1] - 1;
if (igz < 0) igz = 0;
if (igz >= procgrid[2]) igz = procgrid[2] - 1;
return grid2proc[igx][igy][igz];
}
/* ----------------------------------------------------------------------
binary search for value in N-length ascending vec
value may be outside range of vec limits
always return index from 0 to N-1 inclusive
return 0 if value < vec[0]
reutrn N-1 if value >= vec[N-1]
return index = 1 to N-2 if vec[index] <= value < vec[index+1]
------------------------------------------------------------------------- */
int Comm::binary(double value, int n, double *vec)
{
int lo = 0;
int hi = n-1;
if (value < vec[lo]) return lo;
if (value >= vec[hi]) return hi;
// insure vec[lo] <= value < vec[hi] at every iteration
// done when lo,hi are adjacent
int index = (lo+hi)/2;
while (lo < hi-1) {
if (value < vec[index]) hi = index;
else if (value >= vec[index]) lo = index;
index = (lo+hi)/2;
}
return index;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
communicate inbuf around full ring of processors with messtag communicate inbuf around full ring of processors with messtag
nbytes = size of inbuf = n datums * nper bytes nbytes = size of inbuf = n datums * nper bytes

View File

@ -21,9 +21,8 @@ namespace LAMMPS_NS {
class Comm : protected Pointers { class Comm : protected Pointers {
public: public:
int style; // comm pattern: 0 = 6-way stencil, 1 = irregular tiling int style; // comm pattern: 0 = 6-way stencil, 1 = irregular tiling
int layout; // LAYOUT_UNIFORM = logical equal-sized bricks int layout; // LAYOUT_UNIFORM = equal-sized bricks
// LAYOUT_NONUNIFORM = logical bricks, // LAYOUT_NONUNIFORM = logical bricks, but diff sizes via LB
// but different sizes due to LB
// LAYOUT_TILED = general tiling, due to RCB LB // LAYOUT_TILED = general tiling, due to RCB LB
int me,nprocs; // proc info int me,nprocs; // proc info
@ -89,6 +88,15 @@ class Comm : protected Pointers {
virtual void forward_comm_array(int, double **) = 0; virtual void forward_comm_array(int, double **) = 0;
virtual int exchange_variable(int, double *, double *&) = 0; virtual int exchange_variable(int, double *, double *&) = 0;
int binary(double, int, double *);
// map a point to a processor, based on current decomposition
virtual void coord2proc_setup() {}
virtual int coord2proc(double *, int &, int &, int &);
// memory usage
virtual bigint memory_usage() = 0; virtual bigint memory_usage() = 0;
// non-virtual functions common to all Comm styles // non-virtual functions common to all Comm styles

View File

@ -159,30 +159,33 @@ void CommTiled::setup()
// if RCB decomp exists and just changed, gather needed global RCB info // if RCB decomp exists and just changed, gather needed global RCB info
if (rcbnew) { if (layout == LAYOUT_TILED) coord2proc_setup();
if (!rcbinfo)
rcbinfo = (RCBinfo *)
memory->smalloc(nprocs*sizeof(RCBinfo),"comm:rcbinfo");
rcbnew = 0;
RCBinfo rcbone;
memcpy(&rcbone.mysplit[0][0],&mysplit[0][0],6*sizeof(double));
rcbone.cutfrac = rcbcutfrac;
rcbone.dim = rcbcutdim;
MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR,
rcbinfo,sizeof(RCBinfo),MPI_CHAR,world);
}
// set cutoff for comm forward and comm reverse
// check that cutoff < any periodic box length // check that cutoff < any periodic box length
double cut = MAX(neighbor->cutneighmax,cutghostuser); double cut = MAX(neighbor->cutneighmax,cutghostuser);
cutghost[0] = cutghost[1] = cutghost[2] = cut; cutghost[0] = cutghost[1] = cutghost[2] = cut;
if ((periodicity[0] && cut > prd[0]) || if ((periodicity[0] && cut > prd[0]) ||
(periodicity[1] && cut > prd[1]) || (periodicity[1] && cut > prd[1]) ||
(dimension == 3 && periodicity[2] && cut > prd[2])) (dimension == 3 && periodicity[2] && cut > prd[2]))
error->all(FLERR,"Communication cutoff for comm_style tiled " error->all(FLERR,"Communication cutoff for comm_style tiled "
"cannot exceed periodic box length"); "cannot exceed periodic box length");
// if cut = 0.0, set to epsilon to induce nearest neighbor comm
// this is b/c sendproc is used below to infer touching exchange procs
// exchange procs will be empty (leading to lost atoms) if sendproc = 0
// will reset sendproc/etc to 0 after exchange is setup, down below
int cutzero = 0;
if (cut == 0.0) {
cutzero = 1;
cut = MIN(prd[0],prd[1]);
if (dimension == 3) cut = MIN(cut,prd[3]);
cut *= EPSILON*EPSILON;
}
// setup forward/reverse communication // setup forward/reverse communication
// loop over 6 swap directions // loop over 6 swap directions
// determine which procs I will send to and receive from in each swap // determine which procs I will send to and receive from in each swap
@ -416,6 +419,15 @@ void CommTiled::setup()
for (i = 0; i < noverlap; i++) exchproc[idim][i] = overlap[i]; for (i = 0; i < noverlap; i++) exchproc[idim][i] = overlap[i];
} }
// reset sendproc/etc to 0 if cut is really 0.0
if (cutzero) {
for (i = 0; i < nswap; i++) {
nsendproc[i] = nrecvproc[i] =
sendother[i] = recvother[i] = sendself[i] = 0;
}
}
// reallocate MPI Requests and Statuses as needed // reallocate MPI Requests and Statuses as needed
int nmax = 0; int nmax = 0;
@ -999,13 +1011,15 @@ void CommTiled::forward_comm_pair(Pair *pair)
nsize*recvnum[iswap][i], nsize*recvnum[iswap][i],
MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
} }
if (sendother[iswap]) { if (sendother[iswap]) {
for (i = 0; i < nsendproc[iswap]; i++) { for (i = 0; i < nsend; i++) {
n = pair->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i], n = pair->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i],
buf_send,pbc_flag[iswap][i],pbc[iswap][i]); buf_send,pbc_flag[iswap][i],pbc[iswap][i]);
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world);
} }
} }
if (sendself[iswap]) { if (sendself[iswap]) {
pair->pack_forward_comm(sendnum[iswap][nsend],sendlist[iswap][nsend], pair->pack_forward_comm(sendnum[iswap][nsend],sendlist[iswap][nsend],
buf_send,pbc_flag[iswap][nsend], buf_send,pbc_flag[iswap][nsend],
@ -1098,7 +1112,7 @@ void CommTiled::forward_comm_fix(Fix *fix, int size)
MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
} }
if (sendother[iswap]) { if (sendother[iswap]) {
for (i = 0; i < nsendproc[iswap]; i++) { for (i = 0; i < nsend; i++) {
n = fix->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i], n = fix->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i],
buf_send,pbc_flag[iswap][i],pbc[iswap][i]); buf_send,pbc_flag[iswap][i],pbc[iswap][i]);
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world);
@ -1196,7 +1210,7 @@ void CommTiled::forward_comm_compute(Compute *compute)
MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
} }
if (sendother[iswap]) { if (sendother[iswap]) {
for (i = 0; i < nsendproc[iswap]; i++) { for (i = 0; i < nsend; i++) {
n = compute->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i], n = compute->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i],
buf_send,pbc_flag[iswap][i], buf_send,pbc_flag[iswap][i],
pbc[iswap][i]); pbc[iswap][i]);
@ -1292,7 +1306,7 @@ void CommTiled::forward_comm_dump(Dump *dump)
MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
} }
if (sendother[iswap]) { if (sendother[iswap]) {
for (i = 0; i < nsendproc[iswap]; i++) { for (i = 0; i < nsend; i++) {
n = dump->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i], n = dump->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i],
buf_send,pbc_flag[iswap][i], buf_send,pbc_flag[iswap][i],
pbc[iswap][i]); pbc[iswap][i]);
@ -1394,7 +1408,7 @@ void CommTiled::forward_comm_array(int nsize, double **array)
MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
} }
if (sendother[iswap]) { if (sendother[iswap]) {
for (i = 0; i < nsendproc[iswap]; i++) { for (i = 0; i < nsend; i++) {
m = 0; m = 0;
for (iatom = 0; iatom < sendnum[iswap][i]; iatom++) { for (iatom = 0; iatom < sendnum[iswap][i]; iatom++) {
j = sendlist[iswap][i][iatom]; j = sendlist[iswap][i][iatom];
@ -1740,7 +1754,7 @@ int CommTiled::point_drop_tiled(int idim, double *x)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
recursive form recursive point drop thru RCB tree
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int CommTiled::point_drop_tiled_recurse(double *x, int CommTiled::point_drop_tiled_recurse(double *x,
@ -1786,12 +1800,46 @@ int CommTiled::closer_subbox_edge(int idim, double *x)
return 1; return 1;
} }
/* ----------------------------------------------------------------------
if RCB decomp exists and just changed, gather needed global RCB info
------------------------------------------------------------------------- */
void CommTiled::coord2proc_setup()
{
if (!rcbnew) return;
if (!rcbinfo)
rcbinfo = (RCBinfo *)
memory->smalloc(nprocs*sizeof(RCBinfo),"comm:rcbinfo");
rcbnew = 0;
RCBinfo rcbone;
memcpy(&rcbone.mysplit[0][0],&mysplit[0][0],6*sizeof(double));
rcbone.cutfrac = rcbcutfrac;
rcbone.dim = rcbcutdim;
MPI_Allgather(&rcbone,sizeof(RCBinfo),MPI_CHAR,
rcbinfo,sizeof(RCBinfo),MPI_CHAR,world);
}
/* ----------------------------------------------------------------------
determine which proc owns atom with coord x[3] based on current decomp
x will be in box (orthogonal) or lamda coords (triclinic)
if layout = UNIFORM or NONUNIFORM, invoke parent method
if layout = TILED, use point_drop_recurse()
return owning proc ID, ignore igx,igy,igz
------------------------------------------------------------------------- */
int CommTiled::coord2proc(double *x, int &igx, int &igy, int &igz)
{
if (layout != LAYOUT_TILED) return Comm::coord2proc(x,igx,igy,igz);
return point_drop_tiled_recurse(x,0,nprocs-1);
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
realloc the size of the send buffer as needed with BUFFACTOR and bufextra realloc the size of the send buffer as needed with BUFFACTOR and bufextra
if flag = 1, realloc if flag = 1, realloc
if flag = 0, don't need to realloc with copy, just free/malloc if flag = 0, don't need to realloc with copy, just free/malloc
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::grow_send(int n, int flag) void CommTiled::grow_send(int n, int flag)
{ {
maxsend = static_cast<int> (BUFFACTOR * n); maxsend = static_cast<int> (BUFFACTOR * n);

View File

@ -44,6 +44,10 @@ class CommTiled : public Comm {
void forward_comm_array(int, double **); // forward comm of array void forward_comm_array(int, double **); // forward comm of array
int exchange_variable(int, double *, double *&); // exchange on neigh stencil int exchange_variable(int, double *, double *&); // exchange on neigh stencil
void coord2proc_setup();
int coord2proc(double *, int &, int &, int &);
bigint memory_usage(); bigint memory_usage();
private: private:

View File

@ -57,7 +57,7 @@ void DisplaceAtoms::command(int narg, char **arg)
if (igroup == -1) error->all(FLERR,"Could not find displace_atoms group ID"); if (igroup == -1) error->all(FLERR,"Could not find displace_atoms group ID");
int groupbit = group->bitmask[igroup]; int groupbit = group->bitmask[igroup];
int style=-1; int style = -1;
if (strcmp(arg[1],"move") == 0) style = MOVE; if (strcmp(arg[1],"move") == 0) style = MOVE;
else if (strcmp(arg[1],"ramp") == 0) style = RAMP; else if (strcmp(arg[1],"ramp") == 0) style = RAMP;
else if (strcmp(arg[1],"random") == 0) style = RANDOM; else if (strcmp(arg[1],"random") == 0) style = RANDOM;

View File

@ -272,7 +272,7 @@ void FixBalance::rebalance()
// else allow caller's comm->exchange() to do it // else allow caller's comm->exchange() to do it
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (domain->triclinic) domain->x2lamda(atom->nlocal);
if (lbstyle == BISECTION) irregular->migrate_atoms(0,sendproc); if (lbstyle == BISECTION) irregular->migrate_atoms(0,1,sendproc);
else if (irregular->migrate_check()) irregular->migrate_atoms(); else if (irregular->migrate_check()) irregular->migrate_atoms();
if (domain->triclinic) domain->lamda2x(atom->nlocal); if (domain->triclinic) domain->lamda2x(atom->nlocal);

View File

@ -92,12 +92,13 @@ Irregular::~Irregular()
unlike exchange(), allows atoms to have moved arbitrarily long distances unlike exchange(), allows atoms to have moved arbitrarily long distances
sets up irregular plan, invokes it, destroys it sets up irregular plan, invokes it, destroys it
sortflag = flag for sorting order of received messages by proc ID sortflag = flag for sorting order of received messages by proc ID
procassign = non-NULL if already know procs atoms are assigned to (from RCB) preassign = 1 if already know procs that atoms are assigned to via RCB
procassign = list of proc assignments for each owned atom
atoms MUST be remapped to be inside simulation box before this is called atoms MUST be remapped to be inside simulation box before this is called
for triclinic: atoms must be in lamda coords (0-1) before this is called for triclinic: atoms must be in lamda coords (0-1) before this is called
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Irregular::migrate_atoms(int sortflag, int *procassign) void Irregular::migrate_atoms(int sortflag, int preassign, int *procassign)
{ {
// clear global->local map since atoms move to new procs // clear global->local map since atoms move to new procs
// clear old ghosts so map_set() at end will operate only on local atoms // clear old ghosts so map_set() at end will operate only on local atoms
@ -109,7 +110,6 @@ void Irregular::migrate_atoms(int sortflag, int *procassign)
atom->avec->clear_bonus(); atom->avec->clear_bonus();
// subbox bounds for orthogonal or triclinic box // subbox bounds for orthogonal or triclinic box
// other comm/domain data used by coord2proc()
double *sublo,*subhi; double *sublo,*subhi;
if (triclinic == 0) { if (triclinic == 0) {
@ -120,19 +120,14 @@ void Irregular::migrate_atoms(int sortflag, int *procassign)
subhi = domain->subhi_lamda; subhi = domain->subhi_lamda;
} }
layout = comm->layout; // if Comm will be called to assign new atom coords to procs,
xsplit = comm->xsplit; // may need to setup RCB info
ysplit = comm->ysplit;
zsplit = comm->zsplit;
boxlo = domain->boxlo;
prd = domain->prd;
procgrid = comm->procgrid; if (!preassign) comm->coord2proc_setup();
grid2proc = comm->grid2proc;
// loop over atoms, flag any that are not in my sub-box // loop over atoms, flag any that are not in my sub-box
// fill buffer with atoms leaving my box, using < and >= // fill buffer with atoms leaving my box, using < and >=
// assign which proc it belongs to via coord2proc() // assign which proc it belongs to via Comm::coord2proc()
// if coord2proc() returns me, due to round-off // if coord2proc() returns me, due to round-off
// in triclinic x2lamda(), then keep atom and don't send // in triclinic x2lamda(), then keep atom and don't send
// when atom is deleted, fill it in with last atom // when atom is deleted, fill it in with last atom
@ -154,7 +149,7 @@ void Irregular::migrate_atoms(int sortflag, int *procassign)
int nsendatom = 0; int nsendatom = 0;
int i = 0; int i = 0;
if (procassign) { if (preassign) {
while (i < nlocal) { while (i < nlocal) {
if (procassign[i] == me) i++; if (procassign[i] == me) i++;
else { else {
@ -174,7 +169,7 @@ void Irregular::migrate_atoms(int sortflag, int *procassign)
if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] || if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] ||
x[i][1] < sublo[1] || x[i][1] >= subhi[1] || x[i][1] < sublo[1] || x[i][1] >= subhi[1] ||
x[i][2] < sublo[2] || x[i][2] >= subhi[2]) { x[i][2] < sublo[2] || x[i][2] >= subhi[2]) {
mproclist[nsendatom] = coord2proc(x[i],igx,igy,igz); mproclist[nsendatom] = comm->coord2proc(x[i],igx,igy,igz);
if (mproclist[nsendatom] == me) i++; if (mproclist[nsendatom] == me) i++;
else { else {
if (nsend > maxsend) grow_send(nsend,1); if (nsend > maxsend) grow_send(nsend,1);
@ -211,6 +206,7 @@ void Irregular::migrate_atoms(int sortflag, int *procassign)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
check if any atoms need to migrate further than one proc away in any dim check if any atoms need to migrate further than one proc away in any dim
if not, caller can decide to use comm->exchange() instead if not, caller can decide to use comm->exchange() instead
should not be called for layout = TILED
atoms must be remapped to be inside simulation box before this is called atoms must be remapped to be inside simulation box before this is called
for triclinic: atoms must be in lamda coords (0-1) before this is called for triclinic: atoms must be in lamda coords (0-1) before this is called
return 1 if migrate required, 0 if not return 1 if migrate required, 0 if not
@ -224,7 +220,6 @@ int Irregular::migrate_check()
if (comm->layout == LAYOUT_TILED) return 1; if (comm->layout == LAYOUT_TILED) return 1;
// subbox bounds for orthogonal or triclinic box // subbox bounds for orthogonal or triclinic box
// other comm/domain data used by coord2proc()
double *sublo,*subhi; double *sublo,*subhi;
if (triclinic == 0) { if (triclinic == 0) {
@ -235,18 +230,8 @@ int Irregular::migrate_check()
subhi = domain->subhi_lamda; subhi = domain->subhi_lamda;
} }
layout = comm->layout;
xsplit = comm->xsplit;
ysplit = comm->ysplit;
zsplit = comm->zsplit;
boxlo = domain->boxlo;
prd = domain->prd;
procgrid = comm->procgrid;
grid2proc = comm->grid2proc;
// loop over atoms, check for any that are not in my sub-box // loop over atoms, check for any that are not in my sub-box
// assign which proc it belongs to via coord2proc() // assign which proc it belongs to via Comm::coord2proc()
// if logical igx,igy,igz of newproc > one away from myloc, set flag = 1 // if logical igx,igy,igz of newproc > one away from myloc, set flag = 1
// this check needs to observe PBC // this check needs to observe PBC
// cannot check via comm->procneigh since it ignores PBC // cannot check via comm->procneigh since it ignores PBC
@ -264,7 +249,7 @@ int Irregular::migrate_check()
if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] || if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] ||
x[i][1] < sublo[1] || x[i][1] >= subhi[1] || x[i][1] < sublo[1] || x[i][1] >= subhi[1] ||
x[i][2] < sublo[2] || x[i][2] >= subhi[2]) { x[i][2] < sublo[2] || x[i][2] >= subhi[2]) {
coord2proc(x[i],igx,igy,igz); comm->coord2proc(x[i],igx,igy,igz);
glo = myloc[0] - 1; glo = myloc[0] - 1;
ghi = myloc[0] + 1; ghi = myloc[0] + 1;
@ -779,80 +764,6 @@ void Irregular::destroy_data()
delete [] status; delete [] status;
} }
/* ----------------------------------------------------------------------
determine which proc owns atom with coord x[3]
x will be in box (orthogonal) or lamda coords (triclinic)
if layout = UNIFORM, calculate owning proc directly
else layout = NONUNIFORM, iteratively find owning proc via binary search
return owning proc ID via grid2proc
return igx,igy,igz = logical grid loc of owing proc within 3d grid of procs
------------------------------------------------------------------------- */
int Irregular::coord2proc(double *x, int &igx, int &igy, int &igz)
{
if (layout == 0) {
if (triclinic == 0) {
igx = static_cast<int> (procgrid[0] * (x[0]-boxlo[0]) / prd[0]);
igy = static_cast<int> (procgrid[1] * (x[1]-boxlo[1]) / prd[1]);
igz = static_cast<int> (procgrid[2] * (x[2]-boxlo[2]) / prd[2]);
} else {
igx = static_cast<int> (procgrid[0] * x[0]);
igy = static_cast<int> (procgrid[1] * x[1]);
igz = static_cast<int> (procgrid[2] * x[2]);
}
} else {
if (triclinic == 0) {
igx = binary((x[0]-boxlo[0])/prd[0],procgrid[0],xsplit);
igy = binary((x[1]-boxlo[1])/prd[1],procgrid[1],ysplit);
igz = binary((x[2]-boxlo[2])/prd[2],procgrid[2],zsplit);
} else {
igx = binary(x[0],procgrid[0],xsplit);
igy = binary(x[1],procgrid[1],ysplit);
igz = binary(x[2],procgrid[2],zsplit);
}
}
if (igx < 0) igx = 0;
if (igx >= procgrid[0]) igx = procgrid[0] - 1;
if (igy < 0) igy = 0;
if (igy >= procgrid[1]) igy = procgrid[1] - 1;
if (igz < 0) igz = 0;
if (igz >= procgrid[2]) igz = procgrid[2] - 1;
return grid2proc[igx][igy][igz];
}
/* ----------------------------------------------------------------------
binary search for value in N-length ascending vec
value may be outside range of vec limits
always return index from 0 to N-1 inclusive
return 0 if value < vec[0]
reutrn N-1 if value >= vec[N-1]
return index = 1 to N-2 if vec[index] <= value < vec[index+1]
------------------------------------------------------------------------- */
int Irregular::binary(double value, int n, double *vec)
{
int lo = 0;
int hi = n-1;
if (value < vec[lo]) return lo;
if (value >= vec[hi]) return hi;
// insure vec[lo] <= value < vec[hi] at every iteration
// done when lo,hi are adjacent
int index = (lo+hi)/2;
while (lo < hi-1) {
if (value < vec[index]) hi = index;
else if (value >= vec[index]) lo = index;
index = (lo+hi)/2;
}
return index;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA
if flag = 1, realloc if flag = 1, realloc

View File

@ -27,7 +27,8 @@ class Irregular : protected Pointers {
Irregular(class LAMMPS *); Irregular(class LAMMPS *);
~Irregular(); ~Irregular();
void migrate_atoms(int sortflag = 0, int *procassign = NULL); void migrate_atoms(int sortflag = 0, int preassign = 0,
int *procassign = NULL);
int migrate_check(); int migrate_check();
int create_data(int, int *, int sortflag = 0); int create_data(int, int *, int sortflag = 0);
void exchange_data(char *, int, char *); void exchange_data(char *, int, char *);
@ -38,12 +39,6 @@ class Irregular : protected Pointers {
int me,nprocs; int me,nprocs;
int triclinic; int triclinic;
int map_style; int map_style;
int layout;
double *xsplit,*ysplit,*zsplit; // ptrs to comm
int *procgrid; // ptr to comm
int ***grid2proc; // ptr to comm
double *boxlo; // ptr to domain
double *prd; // ptr to domain
int maxsend,maxrecv; // size of buf send/recv in # of doubles int maxsend,maxrecv; // size of buf send/recv in # of doubles
double *buf_send,*buf_recv; // bufs used in migrate_atoms double *buf_send,*buf_recv; // bufs used in migrate_atoms
@ -90,7 +85,6 @@ class Irregular : protected Pointers {
void exchange_atom(double *, int *, double *); void exchange_atom(double *, int *, double *);
void destroy_atom(); void destroy_atom();
int coord2proc(double *, int &, int &, int &);
int binary(double, int, double *); int binary(double, int, double *);
void grow_send(int,int); // reallocate send buffer void grow_send(int,int); // reallocate send buffer