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

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 (noverlap == maxoverlap) {
-          if (overlap[j] == proc) break;
+          maxoverlap += DELTA_PROCS;
-        if (j < noverlap1) continue;
+          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<int> (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];
   }

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]);
+}

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

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) {

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