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

This commit is contained in:
sjplimp
2014-08-08 17:55:44 +00:00
parent 9026afed5c
commit c1e12f9508
22 changed files with 299 additions and 316 deletions

View File

@ -1428,9 +1428,6 @@ void MSM::particle_map()
int flag = 0;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int i = 0; i < nlocal; i++) {
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge

View File

@ -309,9 +309,6 @@ void MSMCG::particle_map()
int flag = 0;
int i;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int j = 0; j < num_charged; j++) {
i = is_charged[j];

View File

@ -1875,10 +1875,6 @@ void PPPM::particle_map()
int nlocal = atom->nlocal;
int flag = 0;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int i = 0; i < nlocal; i++) {
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@ -1897,9 +1893,8 @@ void PPPM::particle_map()
if (nx+nlower < nxlo_out || nx+nupper > nxhi_out ||
ny+nlower < nylo_out || ny+nupper > nyhi_out ||
nz+nlower < nzlo_out || nz+nupper > nzhi_out) {
nz+nlower < nzlo_out || nz+nupper > nzhi_out)
flag = 1;
}
}
if (flag) error->one(FLERR,"Out of range atoms - cannot compute PPPM");

View File

@ -282,9 +282,6 @@ void PPPMCG::particle_map()
double **x = atom->x;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;
for (int j = 0; j < num_charged; j++) {
int i = is_charged[j];

View File

@ -4209,9 +4209,6 @@ void PPPMDisp::particle_map(double delx, double dely, double delz,
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;
for (int i = 0; i < nlocal; i++) {

View File

@ -78,9 +78,6 @@ void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz,
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (type[i] == typeO) {

View File

@ -679,9 +679,6 @@ void PPPMStagger::particle_map()
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;
for (int i = 0; i < nlocal; i++) {

View File

@ -73,9 +73,6 @@ void PPPMTIP4P::particle_map()
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (type[i] == typeO) {

View File

@ -263,7 +263,7 @@ void FixBondBreak::post_integrate()
}
commflag = 1;
comm->forward_comm_fix(this);
comm->forward_comm_fix(this,2);
// break bonds
// if both atoms list each other as winning bond partner
@ -342,7 +342,7 @@ void FixBondBreak::post_integrate()
// 1-2 neighs already reflect broken bonds
commflag = 2;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this);
// create list of broken bonds that influence my owned atoms
// even if between owned-ghost or ghost-ghost atoms

View File

@ -305,7 +305,7 @@ void FixBondCreate::setup(int vflag)
// if newton_bond is set, need to sum bondcount
commflag = 1;
if (newton_bond) comm->reverse_comm_fix(this);
if (newton_bond) comm->reverse_comm_fix(this,1);
}
/* ---------------------------------------------------------------------- */
@ -336,7 +336,7 @@ void FixBondCreate::post_integrate()
// forward comm of bondcount, so ghosts have it
commflag = 1;
comm->forward_comm_fix(this);
comm->forward_comm_fix(this,1);
// resize bond partner list and initialize it
// probability array overlays distsq array
@ -448,7 +448,7 @@ void FixBondCreate::post_integrate()
}
commflag = 2;
comm->forward_comm_fix(this);
comm->forward_comm_fix(this,2);
// create bonds for atoms I own
// only if both atoms list each other as winning bond partner
@ -545,7 +545,7 @@ void FixBondCreate::post_integrate()
// 1-2 neighs already reflect created bonds
commflag = 3;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this);
// create list of broken bonds that influence my owned atoms
// even if between owned-ghost or ghost-ghost atoms
@ -1270,7 +1270,6 @@ void FixBondCreate::unpack_forward_comm(int n, int first, double *buf)
}
} else {
int **nspecial = atom->nspecial;
tagint **special = atom->special;

View File

@ -576,7 +576,7 @@ void FixRigidNHSmall::initial_integrate(int vflag)
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,26);
// accumulate translational and rotational kinetic energies
@ -704,7 +704,7 @@ void FixRigidNHSmall::final_integrate()
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// include Langevin thermostat forces and torques
@ -774,7 +774,7 @@ void FixRigidNHSmall::final_integrate()
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// accumulate translational and rotational kinetic energies

View File

@ -587,7 +587,7 @@ void FixRigidSmall::setup(int vflag)
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// virial setup before call to set_v
@ -603,7 +603,7 @@ void FixRigidSmall::setup(int vflag)
}
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// set velocity/rotation of atoms in rigid bodues
@ -669,7 +669,7 @@ void FixRigidSmall::initial_integrate(int vflag)
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,26);
// set coords/orient and velocity/rotation of atoms in rigid bodies
@ -797,7 +797,7 @@ void FixRigidSmall::final_integrate()
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// include Langevin thermostat forces and torques
@ -839,7 +839,7 @@ void FixRigidSmall::final_integrate()
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
@ -915,7 +915,7 @@ void FixRigidSmall::pre_neighbor()
nghost_body = 0;
commflag = FULL_BODY;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this);
reset_atom2body();
//check(4);
@ -1004,7 +1004,7 @@ int FixRigidSmall::dof(int tgroup)
}
commflag = DOF;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,3);
// nall = count0 = # of point particles in each rigid body
// mall = count1 = # of finite-size particles in each rigid body
@ -1776,7 +1776,7 @@ void FixRigidSmall::setup_bodies_static()
nghost_body = 0;
commflag = FULL_BODY;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this);
reset_atom2body();
// compute mass & center-of-mass of each rigid body
@ -1813,7 +1813,7 @@ void FixRigidSmall::setup_bodies_static()
// reverse communicate xcm, mass of all bodies
commflag = XCM_MASS;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,4);
for (ibody = 0; ibody < nlocal_body; ibody++) {
xcm = body[ibody].xcm;
@ -1931,7 +1931,7 @@ void FixRigidSmall::setup_bodies_static()
// reverse communicate inertia tensor of all bodies
commflag = ITENSOR;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// overwrite Cartesian inertia tensor with file values
@ -1996,7 +1996,7 @@ void FixRigidSmall::setup_bodies_static()
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,26);
// displace = initial atom coords in basis of principal axes
// set displace = 0.0 for atoms not in any rigid body
@ -2131,7 +2131,7 @@ void FixRigidSmall::setup_bodies_static()
// reverse communicate inertia tensor of all bodies
commflag = ITENSOR;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// error check that re-computed momemts of inertia match diagonalized ones
// do not do test for bodies with params read from infile
@ -2278,7 +2278,7 @@ void FixRigidSmall::setup_bodies_dynamic()
// reverse communicate vcm, angmom of all bodies
commflag = VCM_ANGMOM;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// normalize velocity of COM
@ -3149,7 +3149,7 @@ void FixRigidSmall::zero_momentum()
// forward communicate of vcm to all ghost copies
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// set velocity of atoms in rigid bodues
@ -3175,7 +3175,7 @@ void FixRigidSmall::zero_rotation()
// forward communicate of omega to all ghost copies
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// set velocity of atoms in rigid bodues

View File

@ -99,7 +99,7 @@ void FixRigidSmallOMP::initial_integrate(int vflag)
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,26);
// set coords/orient and velocity/rotation of atoms in rigid bodies
@ -188,7 +188,7 @@ void FixRigidSmallOMP::final_integrate()
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_variable_fix(this);
comm->reverse_comm_fix(this,6);
// include Langevin thermostat forces and torques
@ -236,7 +236,7 @@ void FixRigidSmallOMP::final_integrate()
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm_variable_fix(this);
comm->forward_comm_fix(this,10);
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate

View File

@ -268,7 +268,7 @@ void Balance::command(int narg, char **arg)
// debug output of initial state
#ifdef BALANCE_DEBUG
if (me == 0 && fp) dumpout(update->ntimestep,fp);
if (outflag) dumpout(update->ntimestep,fp);
#endif
int niter = 0;
@ -328,7 +328,7 @@ void Balance::command(int narg, char **arg)
// output of final result
if (outflag && me == 0) dumpout(update->ntimestep,fp);
if (outflag) dumpout(update->ntimestep,fp);
// reset proc sub-domains
// for either brick or tiled comm style
@ -475,19 +475,59 @@ int *Balance::bisection(int sortflag)
{
if (!rcb) rcb = new RCB(lmp);
// NOTE: lo/hi args could be simulation box or particle bounding box
// if particle bbox, then mysplit needs to be reset to sim box
// NOTE: triclinic needs to be in lamda coords
// 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;
rcb->compute(dim,atom->nlocal,atom->x,NULL,boxlo,boxhi);
// shrink-wrap simulation box around atoms for input to RCB
// leads to better-shaped sub-boxes when atoms are far from box boundaries
double shrink[6],shrinkall[6];
shrink[0] = boxhi[0]; shrink[1] = boxhi[1]; shrink[2] = boxhi[2];
shrink[3] = boxlo[0]; shrink[4] = boxlo[1]; shrink[5] = boxlo[2];
double **x = atom->x;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
shrink[0] = MIN(shrink[0],x[i][0]);
shrink[1] = MIN(shrink[1],x[i][1]);
shrink[2] = MIN(shrink[2],x[i][2]);
shrink[3] = MAX(shrink[3],x[i][0]);
shrink[4] = MAX(shrink[4],x[i][1]);
shrink[5] = MAX(shrink[5],x[i][2]);
}
shrink[3] = -shrink[3]; shrink[4] = -shrink[4]; shrink[5] = -shrink[5];
MPI_Allreduce(shrink,shrinkall,6,MPI_DOUBLE,MPI_MIN,world);
shrinkall[3] = -shrinkall[3];
shrinkall[4] = -shrinkall[4];
shrinkall[5] = -shrinkall[5];
double *shrinklo = &shrinkall[0];
double *shrinkhi = &shrinkall[3];
// invoke RCB
// then invert() to create list of proc assignements for my atoms
//rcb->compute(dim,atom->nlocal,atom->x,NULL,boxlo,boxhi);
rcb->compute(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi);
rcb->invert(sortflag);
// reset RCB lo/hi bounding box to full simulation box as needed
// NOTE: this logic is specific to orthogonal boxes, not triclinic
double *lo = rcb->lo;
double *hi = rcb->hi;
if (lo[0] == shrinklo[0]) lo[0] = boxlo[0];
if (lo[1] == shrinklo[1]) lo[1] = boxlo[1];
if (lo[2] == shrinklo[2]) lo[2] = boxlo[2];
if (hi[0] == shrinkhi[0]) hi[0] = boxhi[0];
if (hi[1] == shrinkhi[1]) hi[1] = boxhi[1];
if (hi[2] == shrinkhi[2]) hi[2] = boxhi[2];
// store RCB cut, dim, lo/hi box in CommTiled
// cut and lo/hi need to be in fractional form so can
@ -502,18 +542,18 @@ int *Balance::bisection(int sortflag)
comm->rcbcutdim = idim;
double (*mysplit)[2] = comm->mysplit;
mysplit[0][0] = (lo[0] - boxlo[0]) / prd[0];
if (hi[0] == boxhi[0]) mysplit[0][1] = 1.0;
else mysplit[0][1] = (hi[0] - boxlo[0]) / prd[0];
mysplit[0][0] = (rcb->lo[0] - boxlo[0]) / prd[0];
if (rcb->hi[0] == boxhi[0]) mysplit[0][1] = 1.0;
else mysplit[0][1] = (rcb->hi[0] - boxlo[0]) / prd[0];
mysplit[1][0] = (lo[1] - boxlo[1]) / prd[1];
if (hi[1] == boxhi[1]) mysplit[1][1] = 1.0;
else mysplit[1][1] = (hi[1] - boxlo[1]) / prd[1];
mysplit[1][0] = (rcb->lo[1] - boxlo[1]) / prd[1];
if (rcb->hi[1] == boxhi[1]) mysplit[1][1] = 1.0;
else mysplit[1][1] = (rcb->hi[1] - boxlo[1]) / prd[1];
mysplit[2][0] = (rcb->lo[2] - boxlo[2]) / prd[2];
if (rcb->hi[2] == boxhi[2]) mysplit[2][1] = 1.0;
else mysplit[2][1] = (rcb->hi[2] - boxlo[2]) / prd[2];
mysplit[2][0] = (lo[2] - boxlo[2]) / prd[2];
if (hi[2] == boxhi[2]) mysplit[2][1] = 1.0;
else mysplit[2][1] = (hi[2] - boxlo[2]) / prd[2];
// return list of procs to send my atoms to
@ -669,7 +709,7 @@ int Balance::shift()
#ifdef BALANCE_DEBUG
if (me == 0) debug_shift_output(idim,m+1,np,split);
if (me == 0 && fp) dumpout(update->ntimestep,fp);
if (outflag) dumpout(update->ntimestep,fp);
#endif
// stop if no change in splits, b/c all targets are met exactly
@ -873,89 +913,36 @@ int Balance::binary(double value, int n, double *vec)
write xy lines around each proc's sub-domain for 2d
write xyz cubes around each proc's sub-domain for 3d
only called by proc 0
NOTE: only implemented for orthogonal boxes, not triclinic
------------------------------------------------------------------------- */
void Balance::dumpout(bigint tstep, FILE *bfp)
void Balance::dumpout(bigint tstep, FILE *fp)
{
int dimension = domain->dimension;
// write out one square/cube per processor for 2d/3d
// only write once since topology is static
// write out nodal coords
// some will be duplicates
if (firststep) {
firststep = 0;
fprintf(bfp,"ITEM: TIMESTEP\n");
fprintf(bfp,BIGINT_FORMAT "\n",tstep);
if (dimension == 2) fprintf(bfp,"ITEM: NUMBER OF SQUARES\n");
else fprintf(bfp,"ITEM: NUMBER OF CUBES\n");
fprintf(bfp,"%d\n",nprocs);
if (dimension == 2) fprintf(bfp,"ITEM: SQUARES\n");
else fprintf(bfp,"ITEM: CUBES\n");
/*
double *sublo = domain->sublo;
double *subhi = domain->subhi;
int nx = comm->procgrid[0] + 1;
int ny = comm->procgrid[1] + 1;
int nz = comm->procgrid[2] + 1;
if (dimension == 2) {
int m = 0;
for (int j = 0; j < comm->procgrid[1]; j++)
for (int i = 0; i < comm->procgrid[0]; i++) {
int c1 = j*nx + i + 1;
int c2 = c1 + 1;
int c3 = c2 + nx;
int c4 = c3 - 1;
fprintf(bfp,"%d %d %d %d %d %d\n",m+1,m+1,c1,c2,c3,c4);
m++;
}
} else {
int m = 0;
for (int k = 0; k < comm->procgrid[2]; k++)
for (int j = 0; j < comm->procgrid[1]; j++)
for (int i = 0; i < comm->procgrid[0]; i++) {
int c1 = k*ny*nx + j*nx + i + 1;
int c2 = c1 + 1;
int c3 = c2 + nx;
int c4 = c3 - 1;
int c5 = c1 + ny*nx;
int c6 = c2 + ny*nx;
int c7 = c3 + ny*nx;
int c8 = c4 + ny*nx;
fprintf(bfp,"%d %d %d %d %d %d %d %d %d %d\n",
m+1,m+1,c1,c2,c3,c4,c5,c6,c7,c8);
m++;
}
}
}
// write out nodal coords, can be different every call
// scale xsplit,ysplit,zsplit values to full box
// only implmented for orthogonal boxes, not triclinic
int nx = comm->procgrid[0] + 1;
int ny = comm->procgrid[1] + 1;
int nz = comm->procgrid[2] + 1;
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
fprintf(bfp,"ITEM: TIMESTEP\n");
fprintf(bfp,BIGINT_FORMAT "\n",tstep);
fprintf(bfp,"ITEM: NUMBER OF NODES\n");
if (dimension == 2) fprintf(bfp,"%d\n",nx*ny);
else fprintf(bfp,"%d\n",nx*ny*nz);
fprintf(bfp,"ITEM: BOX BOUNDS\n");
fprintf(bfp,"%g %g\n",boxlo[0],boxhi[0]);
fprintf(bfp,"%g %g\n",boxlo[1],boxhi[1]);
fprintf(bfp,"%g %g\n",boxlo[2],boxhi[2]);
fprintf(bfp,"ITEM: NODES\n");
fprintf(fp,"ITEM: TIMESTEP\n");
fprintf(fp,BIGINT_FORMAT "\n",tstep);
fprintf(fp,"ITEM: NUMBER OF NODES\n");
if (dimension == 2) fprintf(fp,"%d\n",4*nprocs);
else fprintf(fp,"%d\n",8*nprocs);
fprintf(fp,"ITEM: BOX BOUNDS\n");
fprintf(fp,"%g %g\n",boxlo[0],boxhi[0]);
fprintf(fp,"%g %g\n",boxlo[1],boxhi[1]);
fprintf(fp,"%g %g\n",boxlo[2],boxhi[2]);
fprintf(fp,"ITEM: NODES\n");
if (dimension == 2) {
int m = 0;
for (int j = 0; j < ny; j++)
for (int i = 0; i < nx; i++) {
fprintf(bfp,"%d %d %g %g %g\n",m+1,1,
fprintf(fp,"%d %d %g %g %g\n",m+1,1,
boxlo[0] + prd[0]*comm->xsplit[i],
boxlo[1] + prd[1]*comm->ysplit[j],
0.0);
@ -966,13 +953,59 @@ void Balance::dumpout(bigint tstep, FILE *bfp)
for (int k = 0; k < nz; k++)
for (int j = 0; j < ny; j++)
for (int i = 0; i < nx; i++) {
fprintf(bfp,"%d %d %g %g %g\n",m+1,1,
fprintf(fp,"%d %d %g %g %g\n",m+1,1,
boxlo[0] + prd[0]*comm->xsplit[i],
boxlo[1] + prd[1]*comm->ysplit[j],
boxlo[2] + prd[2]*comm->zsplit[k]);
m++;
}
}
// write out one square/cube per processor for 2d/3d
fprintf(fp,"ITEM: TIMESTEP\n");
fprintf(fp,BIGINT_FORMAT "\n",tstep);
if (dimension == 2) fprintf(fp,"ITEM: NUMBER OF SQUARES\n");
else fprintf(fp,"ITEM: NUMBER OF CUBES\n");
fprintf(fp,"%d\n",nprocs);
if (dimension == 2) fprintf(fp,"ITEM: SQUARES\n");
else fprintf(fp,"ITEM: CUBES\n");
int nx = comm->procgrid[0] + 1;
int ny = comm->procgrid[1] + 1;
int nz = comm->procgrid[2] + 1;
if (dimension == 2) {
int m = 0;
for (int j = 0; j < comm->procgrid[1]; j++)
for (int i = 0; i < comm->procgrid[0]; i++) {
int c1 = j*nx + i + 1;
int c2 = c1 + 1;
int c3 = c2 + nx;
int c4 = c3 - 1;
fprintf(fp,"%d %d %d %d %d %d\n",m+1,m+1,c1,c2,c3,c4);
m++;
}
} else {
int m = 0;
for (int k = 0; k < comm->procgrid[2]; k++)
for (int j = 0; j < comm->procgrid[1]; j++)
for (int i = 0; i < comm->procgrid[0]; i++) {
int c1 = k*ny*nx + j*nx + i + 1;
int c2 = c1 + 1;
int c3 = c2 + nx;
int c4 = c3 - 1;
int c5 = c1 + ny*nx;
int c6 = c2 + ny*nx;
int c7 = c3 + ny*nx;
int c8 = c4 + ny*nx;
fprintf(fp,"%d %d %d %d %d %d %d %d %d %d\n",
m+1,m+1,c1,c2,c3,c4,c5,c6,c7,c8);
m++;
}
}
*/
}
/* ----------------------------------------------------------------------

View File

@ -161,6 +161,34 @@ void Comm::init()
triclinic = domain->triclinic;
map_style = atom->map_style;
// warn if any proc's sub-box is smaller than neigh skin
// since may lead to lost atoms in exchange()
// really should check every exchange() in case box size is shrinking
// but seems overkill to do that
int flag = 0;
if (!triclinic) {
if (domain->subhi[0] - domain->sublo[0] < neighbor->skin) flag = 1;
if (domain->subhi[1] - domain->sublo[1] < neighbor->skin) flag = 1;
if (domain->dimension == 3)
if (domain->subhi[2] - domain->sublo[2] < neighbor->skin) flag = 1;
} else {
double delta = domain->subhi_lamda[0] - domain->sublo_lamda[0];
if (delta*domain->prd[0] < neighbor->skin) flag = 1;
delta = domain->subhi_lamda[1] - domain->sublo_lamda[1];
if (delta*domain->prd[1] < neighbor->skin) flag = 1;
if (domain->dimension == 3) {
delta = domain->subhi_lamda[2] - domain->sublo_lamda[2];
if (delta*domain->prd[2] < neighbor->skin) flag = 1;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && me == 0)
error->warning(FLERR,"Proc sub-domain size < neighbor skin - "
"could lead to lost atoms");
// comm_only = 1 if only x,f are exchanged in forward/reverse comm
// comm_x_only = 0 if ghost_velocity since velocities are added
@ -181,13 +209,13 @@ void Comm::init()
for (int i = 0; i < modify->nfix; i++)
size_border += modify->fix[i]->comm_border;
// maxexchange = max # of datums/atom in exchange communication
// maxforward = # of datums in largest forward communication
// maxreverse = # of datums in largest reverse communication
// per-atom limits for communication
// maxexchange = max # of datums in exchange comm, set in exchange()
// maxforward = # of datums in largest forward comm
// maxreverse = # of datums in largest reverse comm
// query pair,fix,compute,dump for their requirements
// pair style can force reverse comm even if newton off
maxexchange = BUFMIN + maxexchange_fix;
maxforward = MAX(size_forward,size_border);
maxreverse = size_reverse;

View File

@ -76,10 +76,8 @@ class Comm : protected Pointers {
virtual void forward_comm_pair(class Pair *) = 0;
virtual void reverse_comm_pair(class Pair *) = 0;
virtual void forward_comm_fix(class Fix *) = 0;
virtual void reverse_comm_fix(class Fix *) = 0;
virtual void forward_comm_variable_fix(class Fix *) = 0;
virtual void reverse_comm_variable_fix(class Fix *) = 0;
virtual void forward_comm_fix(class Fix *, int size=0) = 0;
virtual void reverse_comm_fix(class Fix *, int size=0) = 0;
virtual void forward_comm_compute(class Compute *) = 0;
virtual void reverse_comm_compute(class Compute *) = 0;
virtual void forward_comm_dump(class Dump *) = 0;

View File

@ -99,6 +99,11 @@ void CommBrick::init_buffers()
multilo = multihi = NULL;
cutghostmulti = NULL;
// bufextra = max size of one exchanged atom
// = allowed overflow of sendbuf in exchange()
// atomvec, fix reset these 2 maxexchange values if needed
// only necessary if their size > BUFEXTRA
maxexchange = maxexchange_atom + maxexchange_fix;
bufextra = maxexchange + BUFEXTRA;
@ -558,7 +563,7 @@ void CommBrick::reverse_comm()
exchange: move atoms to correct processors
atoms exchanged with all 6 stencil neighbors
send out atoms that have left my box, receive ones entering my box
atoms will be lost if not inside some proc's box
atoms will be lost if not inside a stencil proc's box
can happen if atom moves outside of non-periodic bounary
or if atom moves more than one proc away
this routine called before every reneighboring
@ -586,6 +591,7 @@ void CommBrick::exchange()
atom->avec->clear_bonus();
// insure send buf is large enough for single atom
// bufextra = max size of one atom = allowed overflow of sendbuf
// fixes can change per-atom size requirement on-the-fly
int bufextra_old = bufextra;
@ -713,7 +719,7 @@ void CommBrick::borders()
// check atoms between nfirst and nlast
// for first swaps in a dim, check owned and ghost
// for later swaps in a dim, only check newly arrived ghosts
// store sent atom indices in list for use in future timesteps
// store sent atom indices in sendlist for use in future timesteps
x = atom->x;
if (mode == SINGLE) {
@ -739,8 +745,8 @@ void CommBrick::borders()
else sendflag = 1;
// find send atoms according to SINGLE vs MULTI
// all atoms eligible versus atoms in bordergroup
// only need to limit loop to bordergroup for first sends (ineed < 2)
// all atoms eligible versus only atoms in bordergroup
// can only limit loop to bordergroup for first sends (ineed < 2)
// on these sends, break loop in two: owned (in group) and ghost
if (sendflag) {
@ -859,6 +865,7 @@ void CommBrick::borders()
/* ----------------------------------------------------------------------
forward communication invoked by a Pair
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::forward_comm_pair(Pair *pair)
@ -898,7 +905,7 @@ void CommBrick::forward_comm_pair(Pair *pair)
/* ----------------------------------------------------------------------
reverse communication invoked by a Pair
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::reverse_comm_pair(Pair *pair)
@ -937,17 +944,22 @@ void CommBrick::reverse_comm_pair(Pair *pair)
/* ----------------------------------------------------------------------
forward communication invoked by a Fix
n = constant number of datums per atom
size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */
void CommBrick::forward_comm_fix(Fix *fix)
void CommBrick::forward_comm_fix(Fix *fix, int size)
{
int iswap,n;
int iswap,n,nsize;
double *buf;
MPI_Request request;
MPI_Status status;
int nsize = fix->comm_forward;
if (size) nsize = size;
else nsize = fix->comm_forward;
for (iswap = 0; iswap < nswap; iswap++) {
@ -977,17 +989,22 @@ void CommBrick::forward_comm_fix(Fix *fix)
/* ----------------------------------------------------------------------
reverse communication invoked by a Fix
n = constant number of datums per atom
size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */
void CommBrick::reverse_comm_fix(Fix *fix)
void CommBrick::reverse_comm_fix(Fix *fix, int size)
{
int iswap,n;
int iswap,n,nsize;
double *buf;
MPI_Request request;
MPI_Status status;
int nsize = fix->comm_reverse;
if (size) nsize = size;
else nsize = fix->comm_forward;
for (iswap = nswap-1; iswap >= 0; iswap--) {
@ -1014,84 +1031,9 @@ void CommBrick::reverse_comm_fix(Fix *fix)
}
}
/* ----------------------------------------------------------------------
forward communication invoked by a Fix
n = total datums for all atoms, allows for variable number/atom
------------------------------------------------------------------------- */
void CommBrick::forward_comm_variable_fix(Fix *fix)
{
int iswap,n;
double *buf;
MPI_Request request;
MPI_Status status;
for (iswap = 0; iswap < nswap; iswap++) {
// pack buffer
n = fix->pack_forward_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flag[iswap],pbc[iswap]);
// exchange with another proc
// if self, set recv buffer to send buffer
if (sendproc[iswap] != me) {
if (recvnum[iswap])
MPI_Irecv(buf_recv,maxrecv,MPI_DOUBLE,recvproc[iswap],0,
world,&request);
if (sendnum[iswap])
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
if (recvnum[iswap]) MPI_Wait(&request,&status);
buf = buf_recv;
} else buf = buf_send;
// unpack buffer
fix->unpack_forward_comm(recvnum[iswap],firstrecv[iswap],buf);
}
}
/* ----------------------------------------------------------------------
reverse communication invoked by a Fix
n = total datums for all atoms, allows for variable number/atom
------------------------------------------------------------------------- */
void CommBrick::reverse_comm_variable_fix(Fix *fix)
{
int iswap,n;
double *buf;
MPI_Request request;
MPI_Status status;
for (iswap = nswap-1; iswap >= 0; iswap--) {
// pack buffer
n = fix->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);
// exchange with another proc
// if self, set recv buffer to send buffer
if (sendproc[iswap] != me) {
if (sendnum[iswap])
MPI_Irecv(buf_recv,maxrecv,MPI_DOUBLE,sendproc[iswap],0,
world,&request);
if (recvnum[iswap])
MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world);
if (sendnum[iswap]) MPI_Wait(&request,&status);
buf = buf_recv;
} else buf = buf_send;
// unpack buffer
fix->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
}
}
/* ----------------------------------------------------------------------
forward communication invoked by a Compute
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::forward_comm_compute(Compute *compute)
@ -1131,7 +1073,7 @@ void CommBrick::forward_comm_compute(Compute *compute)
/* ----------------------------------------------------------------------
reverse communication invoked by a Compute
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::reverse_comm_compute(Compute *compute)
@ -1170,7 +1112,7 @@ void CommBrick::reverse_comm_compute(Compute *compute)
/* ----------------------------------------------------------------------
forward communication invoked by a Dump
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::forward_comm_dump(Dump *dump)
@ -1210,7 +1152,7 @@ void CommBrick::forward_comm_dump(Dump *dump)
/* ----------------------------------------------------------------------
reverse communication invoked by a Dump
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommBrick::reverse_comm_dump(Dump *dump)

View File

@ -33,10 +33,10 @@ class CommBrick : public Comm {
virtual void forward_comm_pair(class Pair *); // forward comm from a Pair
virtual void reverse_comm_pair(class Pair *); // reverse comm from a Pair
virtual void forward_comm_fix(class Fix *); // forward comm from a Fix
virtual void reverse_comm_fix(class Fix *); // reverse comm from a Fix
virtual void forward_comm_variable_fix(class Fix *); // variable-size variant
virtual void reverse_comm_variable_fix(class Fix *); // variable-size variant
virtual void forward_comm_fix(class Fix *, int size=0);
// forward comm from a Fix
virtual void reverse_comm_fix(class Fix *, int size=0);
// reverse comm from a Fix
virtual void forward_comm_compute(class Compute *); // forward from a Compute
virtual void reverse_comm_compute(class Compute *); // reverse from a Compute
virtual void forward_comm_dump(class Dump *); // forward comm from a Dump

View File

@ -46,8 +46,6 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
{
error->all(FLERR,"Comm_style tiled is not yet supported");
if (lmp->cuda)
error->all(FLERR,"USER-CUDA package does not yet support comm_style tiled");
if (lmp->kokkos)
@ -62,8 +60,6 @@ CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
CommTiled::CommTiled(LAMMPS *lmp, Comm *oldcomm) : Comm(*oldcomm)
{
error->all(FLERR,"Comm_style tiled is not yet supported");
if (lmp->cuda)
error->all(FLERR,"USER-CUDA package does not yet support comm_style tiled");
if (lmp->kokkos)
@ -92,6 +88,11 @@ CommTiled::~CommTiled()
void CommTiled::init_buffers()
{
// bufextra = max size of one exchanged atom
// = allowed overflow of sendbuf in exchange()
// atomvec, fix reset these 2 maxexchange values if needed
// only necessary if their size > BUFEXTRA
maxexchange = maxexchange_atom + maxexchange_fix;
bufextra = maxexchange + BUFEXTRA;
@ -617,10 +618,9 @@ void CommTiled::reverse_comm()
/* ----------------------------------------------------------------------
exchange: move atoms to correct processors
NOTE: need to re-doc this
atoms exchanged with all 6 stencil neighbors
atoms exchanged with procs that touch sub-box in each of 3 dims
send out atoms that have left my box, receive ones entering my box
atoms will be lost if not inside some proc's box
atoms will be lost if not inside a touching proc's box
can happen if atom moves outside of non-periodic bounary
or if atom moves more than one proc away
this routine called before every reneighboring
@ -645,6 +645,7 @@ void CommTiled::exchange()
atom->avec->clear_bonus();
// insure send buf is large enough for single atom
// bufextra = max size of one atom = allowed overflow of sendbuf
// fixes can change per-atom size requirement on-the-fly
int bufextra_old = bufextra;
@ -691,7 +692,7 @@ void CommTiled::exchange()
/*
// DEBUG:
// test if proc is not in exch list, means will lose atom
// could be that should lose atom
// could be that *should* lose atom
int flag = 0;
for (int k = 0; k < nexchproc[dim]; k++)
if (proc == exchproc[k]) flag = 1;
@ -774,7 +775,7 @@ void CommTiled::exchange()
void CommTiled::borders()
{
int i,m,n,irecv,nlast,nsend,nrecv,ncount,ncountall;
int i,m,n,irecv,nlast,nsend,nrecv,ngroup,ncount,ncountall;
double xlo,xhi,ylo,yhi,zlo,zhi;
double *bbox;
double **x;
@ -791,14 +792,13 @@ void CommTiled::borders()
for (int iswap = 0; iswap < nswap; iswap++) {
// find atoms within rectangles using >= and <
// find atoms within sendboxes using >= and <
// hi test with ">" is important b/c don't want to send an atom
// in lower dim (on boundary) that a proc will recv again in higher dim
// for x-dim swaps, check owned atoms
// for yz-dim swaps, check owned and ghost atoms
// store sent atom indices in list for use in future timesteps
// NOTE: assume SINGLE mode, add back in logic for MULTI mode later
// and for ngroup when bordergroup is set
// store sent atom indices in sendlist for use in future timesteps
// NOTE: assume SINGLE mode, add logic for MULTI mode later
x = atom->x;
if (iswap % 2 == 0) nlast = atom->nlocal + atom->nghost;
@ -811,12 +811,32 @@ void CommTiled::borders()
ncount = 0;
for (i = 0; i < nlast; i++) {
if (x[i][0] >= xlo && x[i][0] < xhi &&
x[i][1] >= ylo && x[i][1] < yhi &&
x[i][2] >= zlo && x[i][2] < zhi) {
if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
sendlist[iswap][m][ncount++] = i;
if (!bordergroup) {
for (i = 0; i < nlast; i++) {
if (x[i][0] >= xlo && x[i][0] < xhi &&
x[i][1] >= ylo && x[i][1] < yhi &&
x[i][2] >= zlo && x[i][2] < zhi) {
if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
sendlist[iswap][m][ncount++] = i;
}
}
} else {
ngroup = atom->nfirst;
for (i = 0; i < ngroup; i++) {
if (x[i][0] >= xlo && x[i][0] < xhi &&
x[i][1] >= ylo && x[i][1] < yhi &&
x[i][2] >= zlo && x[i][2] < zhi) {
if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
sendlist[iswap][m][ncount++] = i;
}
}
for (i = atom->nlocal; i < nlast; i++) {
if (x[i][0] >= xlo && x[i][0] < xhi &&
x[i][1] >= ylo && x[i][1] < yhi &&
x[i][2] >= zlo && x[i][2] < zhi) {
if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
sendlist[iswap][m][ncount++] = i;
}
}
}
@ -961,6 +981,7 @@ void CommTiled::borders()
/* ----------------------------------------------------------------------
forward communication invoked by a Pair
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::forward_comm_pair(Pair *pair)
@ -1007,7 +1028,7 @@ void CommTiled::forward_comm_pair(Pair *pair)
/* ----------------------------------------------------------------------
reverse communication invoked by a Pair
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::reverse_comm_pair(Pair *pair)
@ -1053,15 +1074,20 @@ void CommTiled::reverse_comm_pair(Pair *pair)
/* ----------------------------------------------------------------------
forward communication invoked by a Fix
n = constant number of datums per atom
size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */
void CommTiled::forward_comm_fix(Fix *fix)
void CommTiled::forward_comm_fix(Fix *fix, int size)
{
int i,irecv,n,nsend,nrecv;
int i,irecv,n,nsize,nsend,nrecv;
MPI_Status status;
int nsize = fix->comm_forward;
if (size) nsize = size;
else nsize = fix->comm_forward;
for (int iswap = 0; iswap < nswap; iswap++) {
nsend = nsendproc[iswap] - sendself[iswap];
@ -1100,15 +1126,20 @@ void CommTiled::forward_comm_fix(Fix *fix)
/* ----------------------------------------------------------------------
reverse communication invoked by a Fix
n = constant number of datums per atom
size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */
void CommTiled::reverse_comm_fix(Fix *fix)
void CommTiled::reverse_comm_fix(Fix *fix, int size)
{
int i,irecv,n,nsend,nrecv;
int i,irecv,n,nsize,nsend,nrecv;
MPI_Status status;
int nsize = fix->comm_reverse;
if (size) nsize = size;
else nsize = fix->comm_forward;
for (int iswap = nswap-1; iswap >= 0; iswap--) {
nsend = nsendproc[iswap] - sendself[iswap];
@ -1144,30 +1175,9 @@ void CommTiled::reverse_comm_fix(Fix *fix)
}
}
/* ----------------------------------------------------------------------
forward communication invoked by a Fix
n = total datums for all atoms, allows for variable number/atom
NOTE: complicated b/c don't know # to recv a priori
------------------------------------------------------------------------- */
void CommTiled::forward_comm_variable_fix(Fix *fix)
{
// NOTE: these two forward/reverse methods still need to be updated
}
/* ----------------------------------------------------------------------
reverse communication invoked by a Fix
n = total datums for all atoms, allows for variable number/atom
------------------------------------------------------------------------- */
void CommTiled::reverse_comm_variable_fix(Fix *fix)
{
}
/* ----------------------------------------------------------------------
forward communication invoked by a Compute
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::forward_comm_compute(Compute *compute)
@ -1216,7 +1226,7 @@ void CommTiled::forward_comm_compute(Compute *compute)
/* ----------------------------------------------------------------------
reverse communication invoked by a Compute
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::reverse_comm_compute(Compute *compute)
@ -1263,7 +1273,7 @@ void CommTiled::reverse_comm_compute(Compute *compute)
/* ----------------------------------------------------------------------
forward communication invoked by a Dump
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::forward_comm_dump(Dump *dump)
@ -1311,7 +1321,7 @@ void CommTiled::forward_comm_dump(Dump *dump)
/* ----------------------------------------------------------------------
reverse communication invoked by a Dump
n = constant number of datums per atom
nsize used only to set recv buffer limit
------------------------------------------------------------------------- */
void CommTiled::reverse_comm_dump(Dump *dump)
@ -1444,8 +1454,8 @@ int CommTiled::exchange_variable(int n, double *inbuf, double *&outbuf)
void CommTiled::box_drop_brick(int idim, double *lo, double *hi, int &indexme)
{
// NOTE: this is not triclinic compatible
// NOTE: there error messages are internal - should not occur
// can remove at some point
// NOTE: these error messages are internal sanity checks
// should not occur, can be removed at some point
int index,dir;
if (hi[idim] == sublo[idim]) {

View File

@ -33,10 +33,10 @@ class CommTiled : public Comm {
void forward_comm_pair(class Pair *); // forward comm from a Pair
void reverse_comm_pair(class Pair *); // reverse comm from a Pair
void forward_comm_fix(class Fix *); // forward comm from a Fix
void reverse_comm_fix(class Fix *); // reverse comm from a Fix
void forward_comm_variable_fix(class Fix *); // variable-size variant
void reverse_comm_variable_fix(class Fix *); // variable-size variant
virtual void forward_comm_fix(class Fix *, int size=0);
// forward comm from a Fix
virtual void reverse_comm_fix(class Fix *, int size=0);
// reverse comm from a Fix
void forward_comm_compute(class Compute *); // forward from a Compute
void reverse_comm_compute(class Compute *); // reverse from a Compute
void forward_comm_dump(class Dump *); // forward comm from a Dump

View File

@ -73,12 +73,14 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
// optional args
outflag = 0;
int outarg = 0;
fp = NULL;
while (iarg < narg) {
if (strcmp(arg[iarg],"out") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix balance command");
outflag = 1;
outarg = iarg+1;
iarg += 2;
} else error->all(FLERR,"Illegal fix balance command");
@ -114,7 +116,7 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
// output file
if (outarg && comm->me == 0) {
if (outflag && comm->me == 0) {
fp = fopen(arg[outarg],"w");
if (fp == NULL) error->one(FLERR,"Cannot open fix balance output file");
}
@ -256,9 +258,9 @@ void FixBalance::rebalance()
comm->layout = LAYOUT_TILED;
}
// output of final result
// output of new decomposition
if (fp) balance->dumpout(update->ntimestep,fp);
if (outflag) balance->dumpout(update->ntimestep,fp);
// reset proc sub-domains
@ -269,12 +271,9 @@ void FixBalance::rebalance()
// only needed if migrate_check() says an atom moves to far
// else allow caller's comm->exchange() to do it
//NOTE: change back to migrate_check()
if (domain->triclinic) domain->x2lamda(atom->nlocal);
if (lbstyle == BISECTION) irregular->migrate_atoms(0,sendproc);
//else if (irregular->migrate_check()) irregular->migrate_atoms();
else irregular->migrate_atoms();
else if (irregular->migrate_check()) irregular->migrate_atoms();
if (domain->triclinic) domain->lamda2x(atom->nlocal);
// invoke KSpace setup_grid() to adjust to new proc sub-domains

View File

@ -40,7 +40,7 @@ class FixBalance : public Fix {
double memory_usage();
private:
int nevery,lbstyle,nitermax;
int nevery,lbstyle,nitermax,outflag;
double thresh,stopthresh;
char bstr[3];
FILE *fp;