diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index fff663120b..75cd116649 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -417,24 +417,39 @@ void CommTiled::setup() noverlap = 0; iswap = 2*idim; n = nsendproc[iswap]; - for (i = 0; i < n; i++) { - proc = sendproc[iswap][i]; - if (proc == me) continue; - if ((this->*box_touch)(proc,idim,0)) overlap[noverlap++] = proc; - } - noverlap1 = noverlap; - iswap = 2*idim; - n = nsendproc[iswap]; for (i = 0; i < n; i++) { proc = sendproc[iswap][i]; if (proc == me) continue; if ((this->*box_touch)(proc,idim,0)) { - for (j = 0; j < noverlap1; j++) - if (overlap[j] == proc) break; - if (j < noverlap1) continue; + if (noverlap == maxoverlap) { + maxoverlap += DELTA_PROCS; + memory->grow(overlap,maxoverlap,"comm:overlap"); + } overlap[noverlap++] = proc; } } + noverlap1 = noverlap; + iswap = 2*idim+1; + n = nsendproc[iswap]; + + MPI_Barrier(world); + + for (i = 0; i < n; i++) { + proc = sendproc[iswap][i]; + if (proc == me) continue; + if ((this->*box_touch)(proc,idim,1)) { + for (j = 0; j < noverlap1; j++) + if (overlap[j] == proc) break; + if (j < noverlap1) continue; + if (noverlap == maxoverlap) { + maxoverlap += DELTA_PROCS; + memory->grow(overlap,maxoverlap,"comm:overlap"); + } + overlap[noverlap++] = proc; + } + } + + MPI_Barrier(world); // reallocate esendproc and erecvproc if needed based on noverlap @@ -447,13 +462,14 @@ void CommTiled::setup() } nexchproc[idim] = noverlap; - for (i = 0; i < noverlap; i++) exchproc[iswap][i] = overlap[i]; + for (i = 0; i < noverlap; i++) exchproc[idim][i] = overlap[i]; } // reallocate MPI Requests and Statuses as needed int nmax = 0; - for (iswap = 0; iswap < nswap; iswap++) nmax = MAX(nmax,nprocmax[iswap]); + for (i = 0; i < nswap; i++) nmax = MAX(nmax,nprocmax[i]); + for (i = 0; i < dimension; i++) nmax = MAX(nmax,nexchprocmax[i]); if (nmax > maxreqstat) { maxreqstat = nmax; delete [] requests; @@ -463,11 +479,13 @@ void CommTiled::setup() } // DEBUG - /* + + MPI_Barrier(world); printf("SUBBOX %d: %g %g: %g %g\n",me,sublo[0],sublo[1],subhi[0],subhi[1]); MPI_Barrier(world); + /* for (i = 0; i < nswap; i++) { if (nsendproc[i] == 1) printf("SETUP SEND %d %d: nsend %d self %d sproc0 %d: " @@ -503,8 +521,25 @@ void CommTiled::setup() i,me,nrecvproc[i],sendother[i],recvproc[i][0],recvproc[i][1]); } - MPI_Barrier(world); */ + + for (i = 0; i < dimension; i++) { + if (nexchproc[i] == 1) + printf("SETUP EXCH %d %d: nexch %d exch %d\n", + i,me,nexchproc[i],exchproc[i][0]); + else if (nexchproc[i] == 2) + printf("SETUP EXCH2 %d %d: nexch %d exch %d %d\n", + i,me,nexchproc[i],exchproc[i][0],exchproc[i][1]); + else if (nexchproc[i] == 3) + printf("SETUP EXCH2 %d %d: nexch %d exch %d %d %d\n", + i,me,nexchproc[i],exchproc[i][0],exchproc[i][1],exchproc[i][2]); + else if (nexchproc[i] == 4) + printf("SETUP EXCH2 %d %d: nexch %d exch %d %d %d %d\n", + i,me,nexchproc[i],exchproc[i][0],exchproc[i][1], + exchproc[i][2],exchproc[i][3]); + } + + MPI_Barrier(world); } /* ---------------------------------------------------------------------- @@ -781,13 +816,18 @@ void CommTiled::exchange() // receiver may receive multiple messages, realloc buf_recv if needed nexch = nexchproc[dim]; + if (!nexch) continue; + + //MPI_Barrier(world); + printf("EXCH dim %d me %d nexch %d %d\n",dim,me,nexch,exchproc[0][0]); + //MPI_Barrier(world); for (m = 0; m < nexch; m++) MPI_Irecv(&exchnum[dim][m],1,MPI_INT, exchproc[dim][m],0,world,&requests[m]); for (m = 0; m < nexch; m++) MPI_Send(&nsend,1,MPI_INT,exchproc[dim][m],0,world); - MPI_Waitall(nrecv,requests,statuses); + MPI_Waitall(nexch,requests,statuses); nrecv = 0; for (m = 0; m < nexch; m++) nrecv += exchnum[dim][m]; @@ -801,7 +841,11 @@ void CommTiled::exchange() } for (m = 0; m < nexch; m++) MPI_Send(buf_send,nsend,MPI_DOUBLE,exchproc[dim][m],0,world); - MPI_Waitall(nrecv,requests,statuses); + MPI_Waitall(nexch,requests,statuses); + + //MPI_Barrier(world); + printf("DONE EXCH dim %d me %d nexch %d %d\n",dim,me,nexch,exchproc[0][0]); + //MPI_Barrier(world); // check incoming atoms to see if I own it and they are in my box // if so, add to my list @@ -820,6 +864,11 @@ void CommTiled::exchange() } m += static_cast (buf_recv[m]); } + + //MPI_Barrier(world); + printf("DONE UNPACK dim %d me %d nexch %d %d\n", + dim,me,nexch,exchproc[0][0]); + //MPI_Barrier(world); } if (atom->firstgroupname) atom->first_reorder(); @@ -1891,7 +1940,7 @@ void CommTiled::deallocate_swap(int n) delete [] nexchproc; delete [] nexchprocmax; - for (int i = 0; i < n; i++) { + for (int i = 0; i < n/2; i++) { delete [] exchproc[i]; delete [] exchnum[i]; } diff --git a/src/domain.cpp b/src/domain.cpp index 13e8a190dc..6a743c3023 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1740,3 +1740,44 @@ void Domain::box_corners() corners[7][0] = 1.0; corners[7][1] = 1.0; corners[7][2] = 1.0; lamda2x(corners[7],corners[7]); } + +/* ---------------------------------------------------------------------- + compute 8 corner pts of my triclinic sub-box + 8 are ordered with x changing fastest, then y, finally z + could be more efficient if just coded with xy,yz,xz explicitly +------------------------------------------------------------------------- */ + +void Domain::subbox_corners() +{ + corners[0][0] = sublo_lamda[0]; corners[0][1] = sublo_lamda[1]; + corners[0][2] = sublo_lamda[2]; + lamda2x(corners[0],corners[0]); + + corners[1][0] = subhi_lamda[0]; corners[1][1] = sublo_lamda[1]; + corners[1][2] = sublo_lamda[2]; + lamda2x(corners[1],corners[1]); + + corners[2][0] = sublo_lamda[0]; corners[2][1] = subhi_lamda[1]; + corners[2][2] = sublo_lamda[2]; + lamda2x(corners[2],corners[2]); + + corners[3][0] = subhi_lamda[0]; corners[3][1] = subhi_lamda[1]; + corners[3][2] = sublo_lamda[2]; + lamda2x(corners[3],corners[3]); + + corners[4][0] = sublo_lamda[0]; corners[4][1] = sublo_lamda[1]; + corners[4][2] = subhi_lamda[2]; + lamda2x(corners[4],corners[4]); + + corners[5][0] = subhi_lamda[0]; corners[5][1] = sublo_lamda[1]; + corners[5][2] = subhi_lamda[2]; + lamda2x(corners[5],corners[5]); + + corners[6][0] = sublo_lamda[0]; corners[6][1] = subhi_lamda[1]; + corners[6][2] = subhi_lamda[2]; + lamda2x(corners[6],corners[6]); + + corners[7][0] = subhi_lamda[0]; corners[7][1] = subhi_lamda[1]; + corners[7][2] = subhi_lamda[2]; + lamda2x(corners[7],corners[7]); +} diff --git a/src/domain.h b/src/domain.h index d9ffb721b0..61a2955b98 100644 --- a/src/domain.h +++ b/src/domain.h @@ -128,6 +128,7 @@ class Domain : protected Pointers { void x2lamda(double *, double *, double *, double *); void bbox(double *, double *, double *, double *); void box_corners(); + void subbox_corners(); // minimum image convention check // return 1 if any distance > 1/2 of box size diff --git a/src/dump_image.cpp b/src/dump_image.cpp index 1c2326b474..d5b613c162 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -122,6 +122,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : boxflag = YES; boxdiam = 0.02; axesflag = NO; + subboxflag = NO; // parse optional args @@ -285,6 +286,15 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal dump image command"); iarg += 4; + } else if (strcmp(arg[iarg],"subbox") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); + if (strcmp(arg[iarg+1],"yes") == 0) subboxflag = YES; + else if (strcmp(arg[iarg+1],"no") == 0) subboxflag = NO; + else error->all(FLERR,"Illegal dump image command"); + subboxdiam = force->numeric(FLERR,arg[iarg+2]); + if (subboxdiam < 0.0) error->all(FLERR,"Illegal dump image command"); + iarg += 3; + } else if (strcmp(arg[iarg],"shiny") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); double shiny = force->numeric(FLERR,arg[iarg+1]); @@ -826,6 +836,36 @@ void DumpImage::create_image() } } + // render outline of my sub-box, orthogonal or triclinic + + if (subboxflag) { + double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo); + if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo); + diameter *= boxdiam; + + double *sublo = domain->sublo; + double *subhi = domain->subhi; + + double (*boxcorners)[3]; + double box[8][3]; + if (domain->triclinic == 0) { + box[0][0] = sublo[0]; box[0][1] = sublo[1]; box[0][2] = sublo[2]; + box[1][0] = subhi[0]; box[1][1] = sublo[1]; box[1][2] = sublo[2]; + box[2][0] = sublo[0]; box[2][1] = subhi[1]; box[2][2] = sublo[2]; + box[3][0] = subhi[0]; box[3][1] = subhi[1]; box[3][2] = sublo[2]; + box[4][0] = sublo[0]; box[4][1] = sublo[1]; box[4][2] = subhi[2]; + box[5][0] = subhi[0]; box[5][1] = sublo[1]; box[5][2] = subhi[2]; + box[6][0] = sublo[0]; box[6][1] = subhi[1]; box[6][2] = subhi[2]; + box[7][0] = subhi[0]; box[7][1] = subhi[1]; box[7][2] = subhi[2]; + boxcorners = box; + } else { + domain->subbox_corners(); + boxcorners = domain->corners; + } + + image->draw_box(boxcorners,diameter); + } + // render outline of simulation box, orthogonal or triclinic if (boxflag) { diff --git a/src/dump_image.h b/src/dump_image.h index fdd61e0fba..72d54e5f9c 100644 --- a/src/dump_image.h +++ b/src/dump_image.h @@ -54,6 +54,8 @@ class DumpImage : public DumpCustom { int zoomvar,perspvar; // index to zoom,persp vars int boxflag,axesflag; // 0/1 for draw box and axes double boxdiam,axeslen,axesdiam; // params for drawing box and axes + int subboxflag; + double subboxdiam; int viewflag; // overall view is static or dynamic