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

Merge pull request #2259 from lammps/rcb-tiled-tri

Browse Source 
enable comm_style tiled and load-balancing to work for triclinic boxes
This commit is contained in:
Axel Kohlmeyer
2020-07-29 21:23:20 -04:00
committed by GitHub
parent b00bfc82c8 dd001e80d2
commit df3d04662d
4 changed files with 96 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/comm_style.rst
View File

@ -55,8 +55,7 @@ commands. The decomposition can be changed via the
Restrictions
""""""""""""
Communication style *tiled* cannot be used with *triclinic* simulation
cells.
None
Related commands
""""""""""""""""

30
src/balance.cpp
View File

@ -569,15 +569,24 @@ int *Balance::bisection(int sortflag)
{
if (!rcb) rcb = new RCB(lmp);
// NOTE: this logic is specific to orthogonal boxes, not triclinic
int dim = domain->dimension;
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
int triclinic = domain->triclinic;
double *boxlo,*boxhi,*prd;
if (triclinic == 0) {
boxlo = domain->boxlo;
boxhi = domain->boxhi;
prd = domain->prd;
} else {
boxlo = domain->boxlo_lamda;
boxhi = domain->boxhi_lamda;
prd = domain->prd_lamda;
}
// shrink-wrap simulation box around atoms for input to RCB
// leads to better-shaped sub-boxes when atoms are far from box boundaries
// if triclinic, do this in lamda coords
double shrink[6],shrinkall[6];
@ -586,6 +595,9 @@ int *Balance::bisection(int sortflag)
double **x = atom->x;
int nlocal = atom->nlocal;
if (triclinic) domain->x2lamda(nlocal);
for (int i = 0; i < nlocal; i++) {
shrink[0] = MIN(shrink[0],x[i][0]);
shrink[1] = MIN(shrink[1],x[i][1]);
@ -621,8 +633,9 @@ int *Balance::bisection(int sortflag)
// invoke RCB
// then invert() to create list of proc assignments for my atoms
// if triclinic, RCB operates on lamda coords
// NOTE: (3/2017) can remove undocumented "old" option at some point
// ditto in rcb.cpp
// ditto in rcb.cpp, or make it an option
if (oldrcb) {
if (wtflag) {
@ -636,6 +649,8 @@ int *Balance::bisection(int sortflag)
} else rcb->compute(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi);
}
if (triclinic) domain->lamda2x(nlocal);
rcb->invert(sortflag);
// reset RCB lo/hi bounding box to full simulation box as needed
@ -1270,6 +1285,9 @@ void Balance::debug_shift_output(int idim, int m, int np, double *split)
else if (bdim[idim] == Z) dim = "Z";
fprintf(stderr,"Dimension %s, Iteration %d\n",dim,m);
fprintf(stderr," Count:");
for (i = 0; i <= np; i++) fmt::print(stderr," {}",count[i]);
fprintf(stderr,"\n");
fprintf(stderr," Sum:");
for (i = 0; i <= np; i++) fmt::print(stderr," {}",sum[i]);
fprintf(stderr,"\n");

83
src/comm_tiled.cpp
View File

@ -67,7 +67,7 @@ CommTiled::~CommTiled()
memory->destroy(buf_send);
memory->destroy(buf_recv);
memory->destroy(overlap);
deallocate_swap(nswap);
deallocate_swap(maxswap);
memory->sfree(rcbinfo);
}
@ -85,8 +85,8 @@ void CommTiled::init_buffers()
maxoverlap = 0;
overlap = NULL;
nswap = 2 * domain->dimension;
allocate_swap(nswap);
maxswap = 6;
allocate_swap(maxswap);
rcbinfo = NULL;
}
@ -97,14 +97,14 @@ void CommTiled::init()
{
Comm::init();
nswap = 2 * domain->dimension;
int bufextra_old = bufextra;
init_exchange();
if (bufextra > bufextra_old) grow_send(maxsend+bufextra,2);
// temporary restrictions
if (triclinic)
error->all(FLERR,"Cannot yet use comm_style tiled with triclinic box");
if (mode == Comm::MULTI)
error->all(FLERR,"Cannot yet use comm_style tiled with multi-mode comm");
}
@ -121,13 +121,21 @@ void CommTiled::setup()
// domain properties used in setup method and methods it calls
dimension = domain->dimension;
int *periodicity = domain->periodicity;
if (triclinic == 0) {
prd = domain->prd;
boxlo = domain->boxlo;
boxhi = domain->boxhi;
sublo = domain->sublo;
subhi = domain->subhi;
int *periodicity = domain->periodicity;
} else {
prd = domain->prd_lamda;
boxlo = domain->boxlo_lamda;
boxhi = domain->boxhi_lamda;
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
// set function pointers
@ -155,11 +163,21 @@ void CommTiled::setup()
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
"will be generated. Atoms may get lost.");
cutghost[0] = cutghost[1] = cutghost[2] = cut;
if (triclinic == 0) cutghost[0] = cutghost[1] = cutghost[2] = cut;
else {
double *h_inv = domain->h_inv;
double length0,length1,length2;
length0 = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
cutghost[0] = cut * length0;
length1 = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
cutghost[1] = cut * length1;
length2 = h_inv[2];
cutghost[2] = cut * length2;
}
if ((periodicity[0] && cut > prd[0]) ||
(periodicity[1] && cut > prd[1]) ||
(dimension == 3 && periodicity[2] && cut > prd[2]))
if ((periodicity[0] && cutghost[0] > prd[0]) ||
(periodicity[1] && cutghost[1] > prd[1]) ||
(dimension == 3 && periodicity[2] && cutghost[2] > prd[2]))
error->all(FLERR,"Communication cutoff for comm_style tiled "
"cannot exceed periodic box length");
@ -174,6 +192,7 @@ void CommTiled::setup()
cut = MIN(prd[0],prd[1]);
if (dimension == 3) cut = MIN(cut,prd[2]);
cut *= EPSILON*EPSILON;
cutghost[0] = cutghost[1] = cutghost[2] = cut;
}
// setup forward/reverse communication
@ -200,11 +219,11 @@ void CommTiled::setup()
lo1[0] = sublo[0]; lo1[1] = sublo[1]; lo1[2] = sublo[2];
hi1[0] = subhi[0]; hi1[1] = subhi[1]; hi1[2] = subhi[2];
if (idir == 0) {
lo1[idim] = sublo[idim] - cut;
lo1[idim] = sublo[idim] - cutghost[idim];
hi1[idim] = sublo[idim];
} else {
lo1[idim] = subhi[idim];
hi1[idim] = subhi[idim] + cut;
hi1[idim] = subhi[idim] + cutghost[idim];
}
two = 0;
@ -306,10 +325,16 @@ void CommTiled::setup()
sbox[3] = MIN(oboxhi[0],hi1[0]);
sbox[4] = MIN(oboxhi[1],hi1[1]);
sbox[5] = MIN(oboxhi[2],hi1[2]);
} else {
pbc_flag[iswap][i] = 1;
if (idir == 0) pbc[iswap][i][idim] = 1;
else pbc[iswap][i][idim] = -1;
if (triclinic) {
if (idim == 1) pbc[iswap][i][5] = pbc[iswap][i][idim];
if (idim == 2) pbc[iswap][i][4] = pbc[iswap][i][3] = pbc[iswap][i][idim];
}
sbox[0] = MAX(oboxlo[0],lo2[0]);
sbox[1] = MAX(oboxlo[1],lo2[1]);
sbox[2] = MAX(oboxlo[2],lo2[2]);
@ -320,21 +345,22 @@ void CommTiled::setup()
if (idir == 0) {
sbox[idim] = sublo[idim];
if (i < noverlap1) sbox[3+idim] = MIN(sbox[3+idim]+cut,subhi[idim]);
else sbox[3+idim] = MIN(sbox[3+idim]-prd[idim]+cut,subhi[idim]);
if (i < noverlap1)
sbox[3+idim] = MIN(sbox[3+idim]+cutghost[idim],subhi[idim]);
else sbox[3+idim] = MIN(sbox[3+idim]-prd[idim]+cutghost[idim],subhi[idim]);
} else {
if (i < noverlap1) sbox[idim] = MAX(sbox[idim]-cut,sublo[idim]);
else sbox[idim] = MAX(sbox[idim]+prd[idim]-cut,sublo[idim]);
if (i < noverlap1) sbox[idim] = MAX(sbox[idim]-cutghost[idim],sublo[idim]);
else sbox[idim] = MAX(sbox[idim]+prd[idim]-cutghost[idim],sublo[idim]);
sbox[3+idim] = subhi[idim];
}
if (idim >= 1) {
if (sbox[0] == oboxlo[0]) sbox[0] -= cut;
if (sbox[3] == oboxhi[0]) sbox[3] += cut;
if (sbox[0] == oboxlo[0]) sbox[0] -= cutghost[0];
if (sbox[3] == oboxhi[0]) sbox[3] += cutghost[0];
}
if (idim == 2) {
if (sbox[1] == oboxlo[1]) sbox[1] -= cut;
if (sbox[4] == oboxhi[1]) sbox[4] += cut;
if (sbox[1] == oboxlo[1]) sbox[1] -= cutghost[1];
if (sbox[4] == oboxhi[1]) sbox[4] += cutghost[1];
}
memcpy(sendbox[iswap][i],sbox,6*sizeof(double));
@ -344,6 +370,7 @@ void CommTiled::setup()
}
}
// setup exchange communication = subset of forward/reverse comm procs
// loop over dimensions
// determine which procs I will exchange with in each dimension
@ -653,14 +680,16 @@ void CommTiled::exchange()
// domain properties used in exchange method and methods it calls
// subbox bounds for orthogonal or triclinic
if (triclinic == 0) {
prd = domain->prd;
boxlo = domain->boxlo;
boxhi = domain->boxhi;
if (triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
prd = domain->prd_lamda;
boxlo = domain->boxlo_lamda;
boxhi = domain->boxhi_lamda;
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
@ -687,6 +716,10 @@ void CommTiled::exchange()
if (proc != me) {
buf_send[nsend++] = proc;
nsend += avec->pack_exchange(i,&buf_send[nsend]);
} else {
// DEBUG statment
// error->warning(FLERR,"Losing atom in CommTiled::exchange() send, "
// "likely bad dynamics");
}
avec->copy(nlocal-1,i,1);
nlocal--;
@ -736,6 +769,9 @@ void CommTiled::exchange()
if (value >= lo && value < hi) {
m += avec->unpack_exchange(&buf_recv[m]);
continue;
} else {
// DEBUG statment
// error->warning(FLERR,"Losing atom in CommTiled::exchange() recv");
}
}
m += static_cast<int> (buf_recv[m]);
@ -785,6 +821,7 @@ void CommTiled::borders()
if (iswap % 2 == 0) nlast = atom->nlocal + atom->nghost;
ncountall = 0;
for (m = 0; m < nsendproc[iswap]; m++) {
bbox = sendbox[iswap][m];
xlo = bbox[0]; ylo = bbox[1]; zlo = bbox[2];

1
src/comm_tiled.h
View File

@ -54,6 +54,7 @@ class CommTiled : public Comm {
private:
int nswap; // # of swaps to perform = 2*dim
int maxswap; // largest nswap can be = 6
// forward/reverse comm info, proc lists include self