Merge branch 'charm' of charmgit:charm into charm
authorAaron Becker <akbecker@gmail.com>
Tue, 21 Dec 2010 03:20:31 +0000 (21:20 -0600)
committerAaron Becker <akbecker@gmail.com>
Tue, 21 Dec 2010 03:20:31 +0000 (21:20 -0600)
19 files changed:
examples/bigsim/sdag/jacobi3d/jacobi3d.C
examples/bigsim/sdag/jacobi3d/jacobi3d.ci
examples/charm++/jacobi2d-sdag/README [new file with mode: 0644]
examples/charm++/jacobi2d-sdag/jacobi2d [deleted file]
examples/charm++/jacobi2d-sdag/jacobi2d.C
examples/charm++/jacobi2d-sdag/jacobi2d.ci
examples/charm++/jacobi3d-sdag/Makefile
examples/charm++/jacobi3d-sdag/jacobi3d.C
examples/charm++/jacobi3d-sdag/jacobi3d.ci
src/ck-ldb/GraphBFTLB.C
src/ck-ldb/ckgraph.C
src/ck-ldb/ckgraph.h
src/libs/ck-libs/ParFUM/ParFUM_Adapt.h
src/libs/ck-libs/ParFUM/ParFUM_Mesh_Modify.h
src/libs/ck-libs/ParFUM/ParFUM_internals.h
src/libs/ck-libs/ParFUM/adapt.C
src/libs/ck-libs/ParFUM/adapt_lock.C
src/libs/ck-libs/ParFUM/mesh_modify.C
src/libs/ck-libs/ParFUM/util.C

index 1dd2f7bed390a304b01f3299cc9ed29caf6ebf8d..98d65ba277f43e9e946d41cd82a16e9cede8da29 100644 (file)
@@ -41,8 +41,6 @@
 /*readonly*/ int num_chare_y;
 /*readonly*/ int num_chare_z;
 
-/*readonly*/ int globalBarrier;
-
 static unsigned long next = 1;
 
 int myrand(int numpes) {
@@ -274,8 +272,10 @@ class Jacobi: public CBase_Jacobi {
     // Send ghost faces to the six neighbors
     void begin_iteration(void) {
       if (thisIndex.x == 0 && thisIndex.y == 0 && thisIndex.z == 0) {
-          CkPrintf("Start of iteration %d\n", iterations);
-          BgPrintf("BgPrint> Start of iteration at %f\n");
+        CkPrintf("Start of iteration %d\n", iterations);
+#if CMK_BLUEGENE_CHARM
+        BgPrintf("BgPrint> Start of iteration at %f\n");
+#endif
       }
       iterations++;
 
@@ -298,28 +298,28 @@ class Jacobi: public CBase_Jacobi {
        for(int k=0; k<blockDimZ; ++k) {
          leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
          rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
-      }
+       }
 
       for(int i=0; i<blockDimX; ++i) 
        for(int k=0; k<blockDimZ; ++k) {
          topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
          bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
-      }
+       }
 
       for(int i=0; i<blockDimX; ++i) 
        for(int j=0; j<blockDimY; ++j) {
          frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
          backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
-      }
+       }
 
       // Send my left face
       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
       // Send my right face
       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
-      // Send my top face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
       // Send my bottom face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
+      thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(bottomMsg);
+      // Send my top face
+      thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(topMsg);
       // Send my front face
       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
       // Send my back face
@@ -335,37 +335,37 @@ class Jacobi: public CBase_Jacobi {
          for(int j=0; j<height; ++j) 
            for(int k=0; k<width; ++k) {
              temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
-         }
+           }
          break;
        case RIGHT:
          for(int j=0; j<height; ++j) 
            for(int k=0; k<width; ++k) {
              temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
-         }
+           }
          break;
-       case TOP:
+       case BOTTOM:
          for(int i=0; i<height; ++i) 
            for(int k=0; k<width; ++k) {
              temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
-         }
+           }
          break;
-       case BOTTOM:
+       case TOP:
          for(int i=0; i<height; ++i) 
            for(int k=0; k<width; ++k) {
              temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
-         }
+           }
          break;
        case FRONT:
          for(int i=0; i<height; ++i) 
            for(int j=0; j<width; ++j) {
-             temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
-         }
+             temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
+           }
          break;
        case BACK:
          for(int i=0; i<height; ++i) 
            for(int j=0; j<width; ++j) {
-             temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
-         }
+             temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
+           }
          break;
        default:
           CkAbort("ERROR\n");
index 5ee6132196137868a4ad3e6ef5ae6f9d4b204e02..817b35b8b9b24313b2b41f219c78bb464dc9bf99 100644 (file)
@@ -12,8 +12,6 @@ mainmodule jacobi3d {
   readonly int num_chare_y;
   readonly int num_chare_z;
 
-  readonly int globalBarrier;
-
   message ghostMsg {
     double gh[];
   };
diff --git a/examples/charm++/jacobi2d-sdag/README b/examples/charm++/jacobi2d-sdag/README
new file mode 100644 (file)
index 0000000..345f932
--- /dev/null
@@ -0,0 +1,19 @@
+Jacobi iteration with a 2D Array.
+
+This code uses a 2-D blocked decomposition of a 2-d array with more
+than one element per chare. The code runs for a specified number of
+iterations, using a reduction to start each successive iteration.
+
+A 5-point stencil pattern is used to update the value for each data
+element. For example, the new value for element X is the current
+value of X plus the current values of its left, right, top, and
+bottom neighbors.
+
+     T
+   L X R
+     B
+
+X'  = (X + L + R + T + B) / 5.0
+
+Boundary conditions are obeyed in this example.  Control is expressed
+using structured dagger.
diff --git a/examples/charm++/jacobi2d-sdag/jacobi2d b/examples/charm++/jacobi2d-sdag/jacobi2d
deleted file mode 100755 (executable)
index 5241598..0000000
Binary files a/examples/charm++/jacobi2d-sdag/jacobi2d and /dev/null differ
index fdae2968def0605cbb0ddc82935258a2a8dc9a1d..e175fd966e7446f34d1ec0a87513938bc2e64d54 100644 (file)
@@ -1,7 +1,10 @@
 /** \file jacobi2d.C
  *  Author: Eric Bohm and Abhinav S Bhatele
- *  This is Abhinav's jacobi3d-sdag cut down to 2d by Eric Bohm
- *  Date Created: Dec 7th, 2010
+ *
+ *  This is jacobi3d-sdag cut down to 2d and fixed to be a correct
+ *  implementation of the finite difference method by Eric Bohm.
+ *
+ *  Date Created: Dec 7th, * 2010
  *
  */
 
 /*readonly*/ int num_chare_x;
 /*readonly*/ int num_chare_y;
 
-/*readonly*/ int globalBarrier;
+/*readonly*/ int maxiterations;
 
 static unsigned long next = 1;
 
-int myrand(int numpes) {
-  next = next * 1103515245 + 12345;
-  return((unsigned)(next/65536) % numpes);
-}
-
-// We want to wrap entries around, and because mod operator % 
-// sometimes misbehaves on negative values. -1 maps to the highest value.
-#define wrap_x(a)      (((a)+num_chare_x)%num_chare_x)
-#define wrap_y(a)      (((a)+num_chare_y)%num_chare_y)
 
 
-#define index(a, b)    ( (a)*(blockDimY+2) + (b) )
+#define index(a, b)    ( (b)*(blockDimX+2) + (a) )
 
-#define MAX_ITER               26
-#define WARM_ITER              5
+#define MAX_ITER               100
+#define WARM_ITER              2
 #define LEFT                   1
 #define RIGHT                  2
 #define TOP                    3
 #define BOTTOM                 4
 #define DIVIDEBY5              0.2
+const double THRESHHOLD =  0.001;
 
 double startTime;
 double endTime;
@@ -51,70 +46,76 @@ double endTime;
  *
  */
 class Main : public CBase_Main {
-  public:
-    CProxy_Jacobi array;
-    int iterations;
-
-    Main(CkArgMsg* m) {
-      if ( (m->argc != 3) && (m->argc != 5) ) {
-        CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
-        CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] \n", m->argv[0]);
-        CkAbort("Abort");
-      }
-
-      // set iteration counter to zero
-      iterations = 0;
+public:
+  CProxy_Jacobi array;
+  int iterations;
+  double maxdifference;
+  Main_SDAG_CODE;
+  Main(CkArgMsg* m) {
+    __sdag_init();
+    if ( (m->argc < 3) || (m->argc > 6)) {
+      CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
+      CkPrintf("OR %s [array_size] [block_size] maxiterations\n", m->argv[0]);
+      CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] \n", m->argv[0]);
+      CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] maxiterations\n", m->argv[0]);
+      CkAbort("Abort");
+    }
 
-      // store the main proxy
-      mainProxy = thisProxy;
+    // set iteration counter to zero
+    iterations = 0;
+    // store the main proxy
+    mainProxy = thisProxy;
        
-      if(m->argc == 3) {
-       arrayDimX = arrayDimY = atoi(m->argv[1]);
-        blockDimX = blockDimY  = atoi(m->argv[2]); 
-      }
-      else if (m->argc == 5) {
-        arrayDimX = atoi(m->argv[1]);
-       arrayDimY = atoi(m->argv[2]);
-        blockDimX = atoi(m->argv[3]); 
-       blockDimY = atoi(m->argv[4]); 
-      }
-
-      if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
-        CkAbort("array_size_X % block_size_X != 0!");
-      if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
-        CkAbort("array_size_Y % block_size_Y != 0!");
-
-      num_chare_x = arrayDimX / blockDimX;
-      num_chare_y = arrayDimY / blockDimY;
-
-      // print info
-      CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
-      CkPrintf("Running Jacobi on %d processors with (%d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y);
-      CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
-      CkPrintf("Block Dimensions: %d %d\n", blockDimX, blockDimY);
+    if(m->argc <= 4) {
+      arrayDimX = arrayDimY = atoi(m->argv[1]);
+      blockDimX = blockDimY = atoi(m->argv[2]);
+    }
+    else if (m->argc >= 5) {
+      arrayDimX = atoi(m->argv[1]);
+      arrayDimY = atoi(m->argv[2]);
+      blockDimX = atoi(m->argv[3]); 
+      blockDimY = atoi(m->argv[4]); 
+    }
+    maxiterations=MAX_ITER;
+    if(m->argc==4)
+      maxiterations=atoi(m->argv[3]); 
+    if(m->argc==6)
+      maxiterations=atoi(m->argv[5]); 
       
-      // make proxy and populate array in one call
-      array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y);
+    if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
+      CkAbort("array_size_X % block_size_X != 0!");
+    if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
+      CkAbort("array_size_Y % block_size_Y != 0!");
+
+    num_chare_x = arrayDimX / blockDimX;
+    num_chare_y = arrayDimY / blockDimY;
+
+    // print info
+    CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
+    CkPrintf("Running Jacobi on %d processors with (%d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y);
+    CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
+    CkPrintf("Block Dimensions: %d %d\n", blockDimX, blockDimY);
+    CkPrintf("max iterations %d\n", maxiterations);
+    CkPrintf("Threshhold %.10g\n", THRESHHOLD);
       
-      // initiate computation
-      array.doStep();
-    }
+    // NOTE: boundary conditions must be set based on values
+      
+    // make proxy and populate array in one call
+    array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y);
 
-    // Each worker reports back to here when it completes an iteration
-    void report() {
-      iterations++;
-      if (iterations <= WARM_ITER) {
-       if (iterations == WARM_ITER)
-         startTime = CmiWallTimer();
-       array.doStep();
-      }
-      else {
-       CkPrintf("Completed %d iterations\n", MAX_ITER-1);
-       endTime = CmiWallTimer();
-       CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
-        CkExit();
-      }
-    }
+    // initiate computation
+    thisProxy.run();
+  }
+
+  void done(bool success) {
+      if(success)
+       CkPrintf("Difference %.10g Satisfied Threshhold %.10g in %d Iterations\n", maxdifference,THRESHHOLD,iterations);
+      else
+       CkPrintf("Completed %d Iterations , Difference %lf fails threshhold\n", iterations,maxdifference);
+      endTime = CmiWallTimer();
+      CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(maxiterations-1-WARM_ITER));
+      CkExit();
+  }
 };
 
 /** \class Jacobi
@@ -125,168 +126,245 @@ class Jacobi: public CBase_Jacobi {
   Jacobi_SDAG_CODE
 
   public:
-    int iterations;
-    int imsg;
-
-    double *temperature;
-    double *new_temperature;
-
-    // Constructor, initialize values
-    Jacobi() {
-      __sdag_init();
-      usesAtSync=CmiTrue;
-
-      int i, j;
-      // allocate a two dimensional array
-      temperature = new double[(blockDimX+2) * (blockDimY+2)];
-      new_temperature = new double[(blockDimX+2) * (blockDimY+2)];
-
-      for(i=0; i<blockDimX+2; ++i) {
-       for(j=0; j<blockDimY+2; ++j) {
-           temperature[index(i, j)] = 0.0;
-       } 
+  double *temperature;
+  double *new_temperature;
+  int imsg;
+  int iterations;
+  int numExpected;
+  int istart,ifinish,jstart,jfinish;
+  double maxdifference;
+  bool leftBound, rightBound, topBound, bottomBound;
+  // Constructor, initialize values
+  Jacobi() {
+    __sdag_init();
+    usesAtSync=CmiTrue;
+
+    int i, j;
+    // allocate a two dimensional array
+    temperature = new double[(blockDimX+2) * (blockDimY+2)];
+    new_temperature = new double[(blockDimX+2) * (blockDimY+2)];
+
+    for(i=0; i<blockDimX+2; ++i) {
+      for(j=0; j<blockDimY+2; ++j) {
+       temperature[index(i, j)] = 0.;
+      } 
+    }
+    imsg=0;
+    iterations = 0;
+    numExpected=0;
+    maxdifference=0.;
+    // determine border conditions
+    leftBound=rightBound=topBound=bottomBound=false;
+    istart=jstart=1;
+    ifinish=blockDimX+1;
+    jfinish=blockDimY+1;
+
+    if(thisIndex.x==0)
+      {
+       leftBound=true;
+       istart++;
       }
+    else
+      numExpected++;
 
-      iterations = 0;
-      imsg = 0;
-      constrainBC();
-    }
+    if(thisIndex.x==num_chare_x-1)
+      {
+       rightBound=true;
+       ifinish--;
+      }
+    else
+      numExpected++;
+
+    if(thisIndex.y==0)
+      {
+       topBound=true;
+       jstart++;
+      }
+    else
+      numExpected++;
+
+    if(thisIndex.y==num_chare_y-1)
+      {
+       bottomBound=true;
+       jfinish--;
+      }
+    else
+      numExpected++;
+    constrainBC();
+  }
 
   void pup(PUP::er &p)
   {
     CBase_Jacobi::pup(p);
     __sdag_pup(p);
-    p|iterations;
     p|imsg;
-
+    p|iterations;
+    p|numExpected;
+    p|maxdifference;
+    p|istart; p|ifinish; p|jstart; p|jfinish;
+    p|leftBound; p|rightBound; p|topBound; p|bottomBound;
+    
     size_t size = (blockDimX+2) * (blockDimY+2);
     if (p.isUnpacking()) {
-       temperature = new double[size];
-       new_temperature = new double[size];
-      }
+      temperature = new double[size];
+      new_temperature = new double[size];
+    }
     p(temperature, size);
     p(new_temperature, size);
   }
 
   Jacobi(CkMigrateMessage* m) {__sdag_init();}
 
-    ~Jacobi() { 
-      delete [] temperature; 
-      delete [] new_temperature; 
-    }
-
-    // Send ghost faces to the six neighbors
-    void begin_iteration(void) {
-      AtSync();
-      if (thisIndex.x == 0 && thisIndex.y == 0) {
-          CkPrintf("Start of iteration %d\n", iterations);
-          //BgPrintf("BgPrint> Start of iteration at %f\n");
-      }
-      iterations++;
+  ~Jacobi() { 
+    delete [] temperature; 
+    delete [] new_temperature; 
+  }
 
-      // Copy different faces into messages
-      double *leftGhost =  new double[blockDimY];
-      double *rightGhost =  new double[blockDimY];
-      double *topGhost =  new double[blockDimX];
-      double *bottomGhost =  new double[blockDimX];
+  // Send ghost faces to the six neighbors
+  void begin_iteration(void) {
+    AtSync();
+    if (thisIndex.x == 0 && thisIndex.y == 0) {
+      CkPrintf("Start of iteration %d\n", iterations);
+    }
+    iterations++;
 
-      for(int j=0; j<blockDimY; ++j) {
+    if(!leftBound)
+      {
+       double *leftGhost =  new double[blockDimY];
+       for(int j=0; j<blockDimY; ++j) 
          leftGhost[j] = temperature[index(1, j+1)];
+       thisProxy(thisIndex.x-1, thisIndex.y)
+         .receiveGhosts(iterations, RIGHT, blockDimY, leftGhost);
+       delete [] leftGhost;
+      }
+    if(!rightBound)
+      {
+       double *rightGhost =  new double[blockDimY];
+       for(int j=0; j<blockDimY; ++j) 
          rightGhost[j] = temperature[index(blockDimX, j+1)];
+        thisProxy(thisIndex.x+1, thisIndex.y)
+         .receiveGhosts(iterations, LEFT, blockDimY, rightGhost);
+       delete [] rightGhost;
       }
-
-      for(int i=0; i<blockDimX; ++i) {
+    if(!topBound)
+      {
+       double *topGhost =  new double[blockDimX];
+       for(int i=0; i<blockDimX; ++i) 
          topGhost[i] = temperature[index(i+1, 1)];
-         bottomGhost[i] = temperature[index(i+1, blockDimY)];
-      }
-      // TODO for the inner dimension we can do this in one memcopy
-
-      // Send my left face
-      thisProxy(wrap_x(thisIndex.x-1), thisIndex.y)
-         .receiveGhosts(iterations, RIGHT, blockDimY, leftGhost);
-      // Send my right face
-      thisProxy(wrap_x(thisIndex.x+1), thisIndex.y)
-         .receiveGhosts(iterations, LEFT, blockDimY, rightGhost);
-      // Send my top face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y-1))
+       thisProxy(thisIndex.x, thisIndex.y-1)
          .receiveGhosts(iterations, BOTTOM, blockDimX, topGhost);
-      // Send my bottom face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y+1))
+       delete [] topGhost;
+      }
+    if(!bottomBound)
+      {
+       double *bottomGhost =  new double[blockDimX];
+       for(int i=0; i<blockDimX; ++i) 
+         bottomGhost[i] = temperature[index(i+1, blockDimY)];
+       thisProxy(thisIndex.x, thisIndex.y+1)
          .receiveGhosts(iterations, TOP, blockDimX, bottomGhost);
-    }
+       delete [] bottomGhost;
+      }
+  }
 
-    void processGhosts(int dir, int size, double gh[]) {
-      switch(dir) {
-       case LEFT:
-         for(int j=0; j<size; ++j) {
-             temperature[index(0, j+1)] = gh[j];
-         }
-         break;
-       case RIGHT:
-         for(int j=0; j<size; ++j) {
-             temperature[index(blockDimX+1, j+1)] = gh[j];
-         }
-         break;
-       case TOP:
-         for(int i=0; i<size; ++i) {
-             temperature[index(i+1, 0)] = gh[i];
-         }
-         break;
-       case BOTTOM:
-         for(int i=0; i<size; ++i) {
-             temperature[index(i+1, blockDimY+1)] = gh[i];
-         }
-         break;
-        default:
-          CkAbort("ERROR\n");
+  void processGhosts(int dir, int size, double gh[]) {
+    switch(dir) {
+    case LEFT:
+      for(int j=0; j<size; ++j) {
+       temperature[index(0, j+1)] = gh[j];
+      }
+      break;
+    case RIGHT:
+      for(int j=0; j<size; ++j) {
+       temperature[index(blockDimX+1, j+1)] = gh[j];
+      }
+      break;
+    case TOP:
+      for(int i=0; i<size; ++i) {
+       temperature[index(i+1, 0)] = gh[i];
       }
+      break;
+    case BOTTOM:
+      for(int i=0; i<size; ++i) {
+       temperature[index(i+1, blockDimY+1)] = gh[i];
+      }
+      break;
+    default:
+      CkAbort("ERROR\n");
     }
+  }
 
 
-    void check_and_compute() {
-      compute_kernel();
-
-      // calculate error
-      // not being done right now since we are doing a fixed no. of iterations
+  void check_and_compute(CkCallback cb, int numSteps) {
+    compute_kernel();
 
-      double *tmp;
-      tmp = temperature;
-      temperature = new_temperature;
-      new_temperature = tmp;
+    double *tmp;
+    tmp = temperature;
+    temperature = new_temperature;
+    new_temperature = tmp;
+    constrainBC();
 
-      constrainBC();
+    contribute(sizeof(double), &maxdifference, CkReduction::max_double, cb);
+    if (numSteps > 1)
+      doStep(cb, numSteps-1);
+  }
 
-      if (iterations <= WARM_ITER || iterations >= MAX_ITER)
-       contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
-      else
-       doStep();
-    }
 
-    // Check to see if we have received all neighbor values yet
-    // If all neighbor values have been received, we update our values and proceed
-    void compute_kernel() {
+  // When all neighbor values have been received, 
+  // we update our values and proceed
+  void compute_kernel() {
+    double temperatureIth=0.;
+    double difference=0.;
+    maxdifference=0.;
 #pragma unroll    
-      for(int i=1; i<blockDimX+1; ++i) {
-       for(int j=1; j<blockDimY+1; ++j) {
-           // update my value based on the surrounding values
-           new_temperature[index(i, j)] = (temperature[index(i-1, j)] 
-                                           +  temperature[index(i+1, j)]
-                                           +  temperature[index(i, j-1)]
-                                           +  temperature[index(i, j+1)]
-                                           +  temperature[index(i, j)] ) * DIVIDEBY5;
-       }
+    for(int i=istart; i<ifinish; ++i) {
+      for(int j=jstart; j<jfinish; ++j) {
+       // calculate discrete mean value property 5 pt stencil
+       temperatureIth=(temperature[index(i, j)] 
+                       + temperature[index(i-1, j)] 
+                       +  temperature[index(i+1, j)]
+                       +  temperature[index(i, j-1)]
+                       +  temperature[index(i, j+1)]) * 0.2;
+
+       // update relative error
+       difference=temperatureIth - temperature[index(i, j)];
+       // fix sign without fabs overhead
+       if(difference<0) difference*=-1.0; 
+       maxdifference=(maxdifference>difference) ? maxdifference : difference;
+       new_temperature[index(i, j)] = temperatureIth;
       }
     }
+  }
 
-    // Enforce some boundary conditions
-    void constrainBC() {
-      // Heat left and top  of each chare's block
-      for(int i=1; i<blockDimX+1; ++i)
-         temperature[index(i, 1)] = 255.0;
-      for(int j=1; j<blockDimY+1; ++j)
-         temperature[index(1, j)] = 255.0;
-    }
+  // Enforce some boundary conditions
+  void constrainBC() {
+    if(topBound)
+      for(int i=0; i<blockDimX+2; ++i)
+       temperature[index(i, 1)] = 1.0;
+    if(leftBound)
+      for(int j=0; j<blockDimY+2; ++j)
+       temperature[index(1, j)] = 1.0;
+
+    if(bottomBound)
+      for(int i=0; i<blockDimX+2; ++i)
+       temperature[index(i, blockDimY)] = 0.;
+    if(rightBound)
+      for(int j=0; j<blockDimY+2; ++j)
+       temperature[index(blockDimX, j)] = 0.;
 
+  }
+  // for debugging
+ void dumpMatrix(double *matrix)
+  {
+    CkPrintf("[%d,%d]\n",thisIndex.x, thisIndex.y);
+    for(int i=0; i<blockDimX+2;++i)
+      {
+       for(int j=0; j<blockDimY+2;++j)
+         {
+           CkPrintf("%0.3lf ",matrix[index(i,j)]);
+         }
+       CkPrintf("\n");
+      }
+  }
 };
 
 
index 5e4c7aa872f5dae2b6a0002278df51f6b16a0712..ae9b7fe82de52f74a446e3ff56c8e246616daaa3 100644 (file)
@@ -9,11 +9,39 @@ mainmodule jacobi2d {
   readonly int num_chare_x;
   readonly int num_chare_y;
 
-  readonly int globalBarrier;
+  readonly int maxiterations;
 
   mainchare Main {
     entry Main(CkArgMsg *m);
-    entry void report();
+    entry void run() {
+      while (iterations < WARM_ITER) {
+       atomic {
+         array.doStep(CkCallback(CkIndex_Main::report(NULL), thisProxy), 1);
+       }
+       when report(CkReductionMsg *msg) atomic {
+         iterations++;
+         delete msg;
+       }
+      }
+      atomic "start_timed_portion" {
+       startTime = CmiWallTimer();
+       array.doStep(CkCallback(CkIndex_Main::report(NULL), thisProxy),
+                    maxiterations);
+      }
+      while (iterations < maxiterations) {
+       // Each worker reports back to here when it completes an iteration
+       // Asynchronous check on threshhold satisfaction
+       when report(CkReductionMsg *msg) atomic {
+         iterations++;
+         maxdifference=((double *) msg->getData())[0];
+         delete msg;
+         if (maxdifference - THRESHHOLD < 0)
+           done(true);
+       }
+      }
+      atomic { done(false); }
+    };
+    entry void report(CkReductionMsg *m);
   };
 
   array [2D] Jacobi {
@@ -22,11 +50,11 @@ mainmodule jacobi2d {
     entry void receiveGhosts(int iter, int dir, int size,
                              double ghosts[size]);
 
-    entry void doStep() {
+    entry void doStep(CkCallback cb, int numSteps) {
       atomic "begin_iteration" {
        begin_iteration();
       }
-      for(imsg = 0; imsg < 4; imsg++) {
+      for(imsg = 0; imsg < numExpected; imsg++) {
        // "iterations" keeps track of messages across steps
        when
           receiveGhosts[iterations] (int iter, int dir, int size,
@@ -36,7 +64,7 @@ mainmodule jacobi2d {
           }
       }
       atomic "doWork" {
-       check_and_compute();
+       check_and_compute(cb, numSteps);
       }
     };
   };
index 090067512b634ce63cfd29e71728cf9c06b3b0ee..f98bc642f59ffab9200686f98dd1b2a2bcf08ee7 100644 (file)
@@ -1,4 +1,4 @@
-OPTS   = -O3 -DUSE_TOPOMAP=0
+OPTS   = -O3
 CHARMC = ../../../bin/charmc $(OPTS)
 
 OBJS = jacobi3d.o
index 2320b74d8b610f623706d531399c3cce249520ca..2f2666463e6a5dc6f388aa92c51d4c77719d746f 100644 (file)
@@ -2,8 +2,6 @@
  *  Author: Abhinav S Bhatele
  *  Date Created: June 01st, 2009
  *
- *  This does a topological placement for a 3d jacobi.
- *
  *     
  *           *****************
  *        *               *  *
@@ -26,8 +24,6 @@
 #include "jacobi3d.decl.h"
 #include "TopoManager.h"
 
-// See README for documentation
-
 /*readonly*/ CProxy_Main mainProxy;
 /*readonly*/ int arrayDimX;
 /*readonly*/ int arrayDimY;
@@ -41,8 +37,6 @@
 /*readonly*/ int num_chare_y;
 /*readonly*/ int num_chare_z;
 
-/*readonly*/ int globalBarrier;
-
 static unsigned long next = 1;
 
 int myrand(int numpes) {
@@ -56,8 +50,7 @@ int myrand(int numpes) {
 #define wrap_y(a)      (((a)+num_chare_y)%num_chare_y)
 #define wrap_z(a)      (((a)+num_chare_z)%num_chare_z)
 
-//#define USE_3D_ARRAYS                0
-#define index(a, b, c) ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
+#define index(a,b,c)   ((a)+(b)*(blockDimX+2)+(c)*(blockDimX+2)*(blockDimY+2))
 
 #define MAX_ITER               26
 #define WARM_ITER              5
@@ -124,39 +117,7 @@ class Main : public CBase_Main {
       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
 
       // Create new array of worker chares
-#if USE_TOPOMAP
-      CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z);
-      CkPrintf("Topology Mapping is being done ... \n");
-      CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
-      opts.setMap(map);
-      array = CProxy_Jacobi::ckNew(opts);
-#else
       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
-#endif
-
-      TopoManager tmgr;
-      CkArray *jarr = array.ckLocalBranch();
-      int jmap[num_chare_x][num_chare_y][num_chare_z];
-
-      int hops=0, p;
-      for(int i=0; i<num_chare_x; i++)
-       for(int j=0; j<num_chare_y; j++)
-         for(int k=0; k<num_chare_z; k++) {
-           jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
-         }
-
-      for(int i=0; i<num_chare_x; i++)
-       for(int j=0; j<num_chare_y; j++)
-         for(int k=0; k<num_chare_z; k++) {
-           p = jmap[i][j][k];
-           hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
-           hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
-           hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
-           hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
-           hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
-           hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
-         }
-      CkPrintf("Total Hops: %d\n", hops);
 
       //Start the computation
       array.doStep();
@@ -203,13 +164,10 @@ class Jacobi: public CBase_Jacobi {
       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
 
-      for(i=0; i<blockDimX+2; ++i) {
-       for(j=0; j<blockDimY+2; ++j) {
-         for(k=0; k<blockDimZ+2; ++k) {
+      for(k=0; k<blockDimZ+2; ++k)
+       for(j=0; j<blockDimY+2; ++j)
+         for(i=0; i<blockDimX+2; ++i)
            temperature[index(i, j, k)] = 0.0;
-         }
-       } 
-      }
 
       iterations = 0;
       imsg = 0;
@@ -256,23 +214,23 @@ class Jacobi: public CBase_Jacobi {
       double *frontGhost =  new double[blockDimX*blockDimY];
       double *backGhost =  new double[blockDimX*blockDimY];
 
-      for(int j=0; j<blockDimY; ++j) 
-       for(int k=0; k<blockDimZ; ++k) {
+      for(int k=0; k<blockDimZ; ++k)
+       for(int j=0; j<blockDimY; ++j) {
          leftGhost[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
          rightGhost[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
-      }
+       }
 
-      for(int i=0; i<blockDimX; ++i) 
-       for(int k=0; k<blockDimZ; ++k) {
+      for(int k=0; k<blockDimZ; ++k)
+       for(int i=0; i<blockDimX; ++i) {
          topGhost[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
          bottomGhost[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
-      }
+       }
 
-      for(int i=0; i<blockDimX; ++i) 
-       for(int j=0; j<blockDimY; ++j) {
+      for(int j=0; j<blockDimY; ++j)
+       for(int i=0; i<blockDimX; ++i) {
          frontGhost[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
          backGhost[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
-      }
+       }
 
       // Send my left face
       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z)
@@ -280,12 +238,12 @@ class Jacobi: public CBase_Jacobi {
       // Send my right face
       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z)
          .receiveGhosts(iterations, LEFT, blockDimY, blockDimZ, rightGhost);
-      // Send my top face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z)
-         .receiveGhosts(iterations, BOTTOM, blockDimX, blockDimZ, topGhost);
       // Send my bottom face
-      thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z)
+      thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z)
          .receiveGhosts(iterations, TOP, blockDimX, blockDimZ, bottomGhost);
+      // Send my top face
+      thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z)
+         .receiveGhosts(iterations, BOTTOM, blockDimX, blockDimZ, topGhost);
       // Send my front face
       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1))
          .receiveGhosts(iterations, BACK, blockDimX, blockDimY, frontGhost);
@@ -297,40 +255,40 @@ class Jacobi: public CBase_Jacobi {
     void processGhosts(int dir, int height, int width, double gh[]) {
       switch(dir) {
        case LEFT:
-         for(int j=0; j<height; ++j) 
-           for(int k=0; k<width; ++k) {
+         for(int k=0; k<width; ++k)
+           for(int j=0; j<height; ++j) {
              temperature[index(0, j+1, k+1)] = gh[k*height+j];
-         }
+           }
          break;
        case RIGHT:
-         for(int j=0; j<height; ++j) 
-           for(int k=0; k<width; ++k) {
+         for(int k=0; k<width; ++k)
+           for(int j=0; j<height; ++j) {
              temperature[index(blockDimX+1, j+1, k+1)] = gh[k*height+j];
-         }
+           }
          break;
-       case TOP:
-         for(int i=0; i<height; ++i) 
-           for(int k=0; k<width; ++k) {
+       case BOTTOM:
+         for(int k=0; k<width; ++k)
+           for(int i=0; i<height; ++i) {
              temperature[index(i+1, 0, k+1)] = gh[k*height+i];
-         }
+           }
          break;
-       case BOTTOM:
-         for(int i=0; i<height; ++i) 
-           for(int k=0; k<width; ++k) {
+       case TOP:
+         for(int k=0; k<width; ++k)
+           for(int i=0; i<height; ++i) {
              temperature[index(i+1, blockDimY+1, k+1)] = gh[k*height+i];
-         }
+           }
          break;
        case FRONT:
-         for(int i=0; i<height; ++i) 
-           for(int j=0; j<width; ++j) {
-             temperature[index(i+1, j+1, blockDimZ+1)] = gh[j*height+i];
-         }
+         for(int j=0; j<width; ++j)
+           for(int i=0; i<height; ++i) {
+             temperature[index(i+1, j+1, 0)] = gh[j*height+i];
+           }
          break;
        case BACK:
-         for(int i=0; i<height; ++i) 
-           for(int j=0; j<width; ++j) {
-             temperature[index(i+1, j+1, 0)] = gh[j*height+i];
-         }
+         for(int j=0; j<width; ++j)
+           for(int i=0; i<height; ++i) {
+             temperature[index(i+1, j+1, blockDimZ+1)] = gh[j*height+i];
+           }
          break;
        default:
           CkAbort("ERROR\n");
@@ -361,9 +319,9 @@ class Jacobi: public CBase_Jacobi {
     // If all neighbor values have been received, we update our values and proceed
     void compute_kernel() {
 #pragma unroll    
-      for(int i=1; i<blockDimX+1; ++i) {
-       for(int j=1; j<blockDimY+1; ++j) {
-         for(int k=1; k<blockDimZ+1; ++k) {
+      for(int k=1; k<blockDimZ+1; ++k)
+       for(int j=1; j<blockDimY+1; ++j)
+         for(int i=1; i<blockDimX+1; ++i) {
            // update my value based on the surrounding values
            new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
                                            +  temperature[index(i+1, j, k)]
@@ -372,119 +330,23 @@ class Jacobi: public CBase_Jacobi {
                                            +  temperature[index(i, j, k-1)]
                                            +  temperature[index(i, j, k+1)]
                                            +  temperature[index(i, j, k)] ) * DIVIDEBY7;
-         }
-       }
-      }
+         } // end for
     }
 
     // Enforce some boundary conditions
     void constrainBC() {
       // Heat left, top and front faces of each chare's block
-      for(int i=1; i<blockDimX+1; ++i)
-       for(int k=1; k<blockDimZ+1; ++k)
+      for(int k=1; k<blockDimZ+1; ++k)
+       for(int i=1; i<blockDimX+1; ++i)
          temperature[index(i, 1, k)] = 255.0;
-      for(int j=1; j<blockDimY+1; ++j)
-       for(int k=1; k<blockDimZ+1; ++k)
-         temperature[index(1, j, k)] = 255.0;
-      for(int i=1; i<blockDimX+1; ++i)
+      for(int k=1; k<blockDimZ+1; ++k)
        for(int j=1; j<blockDimY+1; ++j)
+         temperature[index(1, j, k)] = 255.0;
+      for(int j=1; j<blockDimY+1; ++j)
+       for(int i=1; i<blockDimX+1; ++i)
          temperature[index(i, j, 1)] = 255.0;
     }
 
 };
 
-/** \class JacobiMap
- *
- */
-
-class JacobiMap : public CkArrayMap {
-  public:
-    int X, Y, Z;
-    int *mapping;
-
-    JacobiMap(int x, int y, int z) {
-      X = x; Y = y; Z = z;
-      mapping = new int[X*Y*Z];
-
-      // we are assuming that the no. of chares in each dimension is a 
-      // multiple of the torus dimension
-
-      TopoManager tmgr;
-      int dimNX, dimNY, dimNZ, dimNT;
-
-      dimNX = tmgr.getDimNX();
-      dimNY = tmgr.getDimNY();
-      dimNZ = tmgr.getDimNZ();
-      dimNT = tmgr.getDimNT();
-
-      // we are assuming that the no. of chares in each dimension is a 
-      // multiple of the torus dimension
-      int numCharesPerPe = X*Y*Z/CkNumPes();
-
-      int numCharesPerPeX = X / dimNX;
-      int numCharesPerPeY = Y / dimNY;
-      int numCharesPerPeZ = Z / dimNZ;
-      int pe = 0, pes = CkNumPes();
-
-#if USE_BLOCK_RNDMAP
-      int used[pes];
-      for(int i=0; i<pes; i++)
-       used[i] = 0;
-#endif
-
-      if(dimNT < 2) {  // one core per node
-       if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
-       for(int i=0; i<dimNX; i++)
-         for(int j=0; j<dimNY; j++)
-           for(int k=0; k<dimNZ; k++)
-           {
-#if USE_BLOCK_RNDMAP
-             pe = myrand(pes); 
-             while(used[pe]!=0) {
-               pe = myrand(pes); 
-             }
-             used[pe] = 1;
-#endif
-
-             for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
-               for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
-                 for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
-#if USE_TOPOMAP
-                   mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
-#elif USE_BLOCK_RNDMAP
-                   mapping[ci*Y*Z + cj*Z + ck] = pe;
-#endif
-                 }
-           }
-      } else {         // multiple cores per node
-       // In this case, we split the chares in the X dimension among the
-       // cores on the same node. The strange thing I figured out is that
-       // doing this in the Z dimension is not as good.
-       numCharesPerPeX /= dimNT;
-       if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
-
-       for(int i=0; i<dimNX; i++)
-         for(int j=0; j<dimNY; j++)
-           for(int k=0; k<dimNZ; k++)
-             for(int l=0; l<dimNT; l++)
-               for(int ci=(dimNT*i+l)*numCharesPerPeX; ci<(dimNT*i+l+1)*numCharesPerPeX; ci++)
-                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
-                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
-                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
-                   }
-      } // end of if
-
-      if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
-    }
-
-    ~JacobiMap() { 
-      delete [] mapping;
-    }
-
-    int procNum(int, const CkArrayIndex &idx) {
-      int *index = (int *)idx.data();
-      return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
-    }
-};
-
 #include "jacobi3d.def.h"
index d42422ef406bd360b8fd3d0bd423a9c113151d83..e744b1e544fdca0286a3c6d9782f5df704219004 100644 (file)
@@ -12,8 +12,6 @@ mainmodule jacobi3d {
   readonly int num_chare_y;
   readonly int num_chare_z;
 
-  readonly int globalBarrier;
-
   mainchare Main {
     entry Main(CkArgMsg *m);
     entry void report();
@@ -44,8 +42,4 @@ mainmodule jacobi3d {
     };
   };
 
-  group JacobiMap : CkArrayMap {
-    entry JacobiMap(int x, int y, int z);
-  };
-
 };
index b9d1b6cf04cb61c1a2c96ae3dcb4f6176cd7abea..824212b7138a7035dd134a7382bac1807b225270 100644 (file)
@@ -60,10 +60,10 @@ void GraphBFTLB::work(LDStats *stats) {
     start = vertexq.front();
     vertexq.pop();
 
-    for(i = 0; i < ogr->vertices[start].edgeList.size(); i++) {
+    for(i = 0; i < ogr->vertices[start].sendToList.size(); i++) {
       // look at all neighbors of a node in the queue and map them while
       // inserting them in the queue (so we can look at their neighbors next)
-      nbr = ogr->vertices[start].edgeList[i].getNeighborId();
+      nbr = ogr->vertices[start].sendToList[i].getNeighborId();
       if(ogr->vertices[nbr].getNewPe() == -1) {
        vertexq.push(nbr);
 
index 28029f3c206c5941167cc61ba222fcaa85aec8f8..c649759d33c9a98901886ef2e6e996c2fb0a5ac2 100644 (file)
@@ -64,7 +64,8 @@ ObjGraph::ObjGraph(BaseLB::LDStats *stats) {
       from = stats->getHash(commData.sender);
       to = stats->getHash(commData.receiver.get_destObj());
 
-      vertices[from].edgeList.push_back(Edge(to, commData.messages, commData.bytes));
+      vertices[from].sendToList.push_back(Edge(to, commData.messages, commData.bytes));
+      vertices[to].recvFromList.push_back(Edge(from, commData.messages, commData.bytes));
     }
   } // end for
 }
index fbc3d044fd0eb4ff2c99f766c38cc12470c1c160..9eb6bf61e7d151c03810be51ea7658a203a12f49 100644 (file)
@@ -79,8 +79,9 @@ class Vertex {
     inline void setNewPe(int _newpe) { newPe = _newpe; }
     inline bool isMigratable() { return migratable; }
 
-    // undirected edges from this vertex to other vertices
-    std::vector<Edge> edgeList;
+    // list of vertices this vertex sends messages to and receives from
+    std::vector<Edge> sendToList;
+    std::vector<Edge> recvFromList;
 
   private:
     int id;            // index in the LDStats array
index 57ac23c640cba1a60e3a9da1d3a06535004fc5dd..2e536a72d2bc6ed840ae4774e141f72835bcf71b 100644 (file)
@@ -160,9 +160,11 @@ class FEM_AdaptL : public FEM_Adapt {
                                    int e2_n3, int n3, int n4);
 
   /// Acquire an element in our ghost layer, turning it into a local element
-  int eatIntoElement(int e);
+  int eatIntoElement(int e, bool aggressive_node_removal=false);
   /// Test the adaptivity system to see if any nodes are locked
   void residualLockTest();
+  /// Release all currently held locks on this partition
+  void unlockAll();
   /// Test the mesh for corruption in connectivity/adjacency
   void structureTest();
 };
index c735ae3b42c122fae190cea5bbac591e8f8fc874..4c4e7eeb451dd7c27660a9be10b94364d0b68e40 100644 (file)
@@ -34,7 +34,7 @@ void FEM_remove_node(FEM_Mesh *m, int node);
 ///Add an element on this mesh with this connectivity
 int FEM_add_element(FEM_Mesh *m, int* conn, int conn_size, int elem_type=0, int chunkNo=-1);
 ///Remove an element on this mesh
-int FEM_remove_element(FEM_Mesh *m, int element, int elem_type=0, int permanent=-1);
+int FEM_remove_element(FEM_Mesh *m, int element, int elem_type=0, int permanent=-1, bool aggressive_node_removal=false);
 ///Purge the element from this mesh (invalidate entry)
 int FEM_purge_element(FEM_Mesh *m, int element, int elem_type=0);
 
@@ -287,6 +287,7 @@ class removeElemMsg : public CMessage_removeElemMsg {
   int elementid;
   int elemtype;
   int permanent;
+  bool aggressive_node_removal;
 };
 
 ///A message to verify if the IDXL entries for a node/element on one chunk is consistent with another chunk
index 0e2a6f4820d9cfcbac2e6390d7c1e2ecb86188ba..e41a11519d999930a08570e1a3292630b2b42bc0 100644 (file)
@@ -2842,7 +2842,7 @@ class FEM_MUtil {
   ///Add the element with this conn (indices, typeOfIndices) on this chunk (called from remote chunk)
   void addGhostElementRemote(FEM_Mesh *m, int chk, int elemType, int *indices, int *typeOfIndex, int connSize);
   ///Remove this element on this chunk (called from remote chunk)
-  void removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent);
+  void removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent, bool aggressive_node_removal=false);
   ///Remove this ghost element on this chunk, also delete some ghost nodes and elements
   void removeGhostElementRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int numGhostIndex, int *ghostIndices, int numGhostRNIndex, int *ghostRNIndices, int numGhostREIndex, int *ghostREIndices, int numSharedIndex, int *sharedIndices);
 
@@ -2851,7 +2851,7 @@ class FEM_MUtil {
   ///Add this node to the shared Idxl list (called from remote chunk)
   void addToSharedList(FEM_Mesh *m, int fromChk, int sharedIdx);
   ///Acquire the element specified by this ghost index
-  int eatIntoElement(int localIdx);
+  int eatIntoElement(int localIdx, bool aggressive_node_removal=false);
 
   ///Does this chunk know about this node index (on any idxl list with 'chk')
   bool knowsAbtNode(int chk, int nodeId);
@@ -2900,6 +2900,8 @@ class FEM_MUtil {
   void verifyIdxlListRemote(FEM_Mesh *m, int fromChk, int fsize, int type);
   ///Test that there are no remaining acquired locks on this mesh
   int residualLockTest(FEM_Mesh *m);
+  /// Release all currently held locks on this partition
+  void unlockAll(FEM_Mesh* m);
 };
 
 
index 9eead87ccc91f6d1b1c2d6462889865953c43b98..04523f15fcbb1a0e6b1e9adbcd646eeafcff7b7b 100644 (file)
@@ -837,9 +837,9 @@ int FEM_Adapt::vertex_split_help(int n, int n1, int n2, int e1, int e3)
 // ======================  END vertex_split ===================
 
 
-int FEM_AdaptL::eatIntoElement(int e)
+int FEM_AdaptL::eatIntoElement(int e, bool aggressive_node_removal)
 {
-    return theMod->fmUtil->eatIntoElement(e);
+    return theMod->fmUtil->eatIntoElement(e, aggressive_node_removal);
 }
 
 
@@ -849,6 +849,11 @@ void FEM_AdaptL::residualLockTest()
 }
 
 
+void FEM_AdaptL::unlockAll() {
+    theMod->fmUtil->unlockAll(theMod->fmMesh);
+}
+
+
 void FEM_AdaptL::structureTest()
 {
     theMod->fmUtil->StructureTest(theMod->fmMesh);
index 5c52162fefa2a32985e682e81f37f776654217f1..2147bcc0b624bb09cd48ea81ab315b822e221caa 100644 (file)
@@ -6,8 +6,7 @@
 #include "ParFUM.h"
 #include "ParFUM_internals.h"
 
-#define DEBUG_LOCKS
-
+//#define DEBUG_LOCKS
 //#define DEBUG_1
 #define ERVAL -1000000000  //might cause a problem if there are 100million nodes
 #define ERVAL1 -1000000001
@@ -17,6 +16,9 @@
  * nodes which specifies if locks for each nodes has been acquired or not
  */
 int FEM_AdaptL::lockNodes(int *gotlocks, int *lockrnodes, int numRNodes, int *lockwnodes, int numWNodes) {
+#ifdef CPSD
+    return true;
+#endif
   bool donelocks = false;
   int numNodes = numRNodes + numWNodes;
   for(int i=0; i<numNodes; i++) gotlocks[i] = 0;
@@ -98,6 +100,9 @@ int FEM_AdaptL::lockNodes(int *gotlocks, int *lockrnodes, int numRNodes, int *lo
 /** Same as above, instead it unlocks the nodes
  */
 int FEM_AdaptL::unlockNodes(int *gotlocks, int *lockrnodes, int numRNodes, int *lockwnodes, int numWNodes) {
+#ifdef CPSD
+    return true;
+#endif
   bool donelocks = false;
   int numNodes = numRNodes + numWNodes;
   int *ungetlocks = (int*)malloc(numNodes*sizeof(int));
@@ -193,7 +198,7 @@ int FEM_AdaptL::edge_flip(int n1, int n2) {
       numtries++;
       if(numtries>=13+(7*theMod->idx)%43) {
        if(!warned) {
-         //CkPrintf("[%d]Warning: Possibly a livelock in edge_flip %d & %d, supporting %d, %d\n",theMod->idx,n1,n2,n3,n4);
+         CkPrintf("[%d]Warning: Possibly a livelock in edge_flip %d & %d, supporting %d, %d\n",theMod->idx,n1,n2,n3,n4);
          warned = true;
        }
        CthYield();
@@ -242,17 +247,14 @@ int FEM_AdaptL::edge_bisect(int n1, int n2) {
   for(int i=0; i<numNodes; i++) {
     gotlocks[i] = -1;
   }
-  while(!done) {
 #ifdef CPSD
-      int gotlock = 1;
+  isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
 #else
+  while(!done) {
     int gotlock = lockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
-#endif
     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
     if(isEdge == -1) {
-#ifndef CPSD
       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
-#endif
 #ifdef DEBUG_1
       CkPrintf("[%d]Warning: Bisect %d->%d not done as it is no longer a valid edge\n",theMod->idx,n1,n2);
 #endif
@@ -262,9 +264,7 @@ int FEM_AdaptL::edge_bisect(int n1, int n2) {
       done = true;
     }
     else {
-#ifndef CPSD
       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
-#endif
       locknodes[2] = n3;
       locknodes[3] = n4;
       numtries++;
@@ -280,11 +280,11 @@ int FEM_AdaptL::edge_bisect(int n1, int n2) {
       CthYield();
     }
   }
+#endif
   int ret = edge_bisect_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, e2_n1, e2_n2, e2_n3, n3, n4);
 #ifndef CPSD
   unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
 #endif
-
   return ret;
 }
 
@@ -368,7 +368,7 @@ int FEM_AdaptL::vertex_remove(int n1, int n2) {
       numtries++;
       if(numtries>=13+(7*theMod->idx)%43) {
        if(!warned) {
-         //CkPrintf("[%d]Warning: Possibly a livelock in vertex_remove %d & %d, supporting %d, %d and %d\n",theMod->idx,n1,n2,n3,n4,n5);
+         CkPrintf("[%d]Warning: Possibly a livelock in vertex_remove %d & %d, supporting %d, %d and %d\n",theMod->idx,n1,n2,n3,n4,n5);
          warned = true;
        }
        numtries = 0;
@@ -622,7 +622,7 @@ int FEM_AdaptL::edge_contraction(int n1, int n2) {
        }
       }
       if(numtries>=50) {
-       //CkPrintf("Possibly a livelock in cloud nodes edge_contract\n");
+       CkPrintf("Possibly a livelock in cloud nodes edge_contract\n");
        //it is ok to skip an edge_contract, if the lock is too difficult to get
        isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
        if(isEdge!=-1) {
@@ -646,7 +646,7 @@ int FEM_AdaptL::edge_contraction(int n1, int n2) {
       numtries++;
       if(numtries>=50) {
        if(!warned) {
-         //CkPrintf("[%d]Warning: Possibly a livelock in edge_contract %d & %d, supporting %d, %d. Avoiding this contract operation.\n",theMod->idx,n1,n2,n3,n4);
+         CkPrintf("[%d]Warning: Possibly a livelock in edge_contract %d & %d, supporting %d, %d. Avoiding this contract operation.\n",theMod->idx,n1,n2,n3,n4);
          warned = true;
        }
         //it is ok to skip an edge_contract, if the lock is too difficult to get
@@ -1199,7 +1199,7 @@ int FEM_AdaptL::edge_contraction_help(int *e1P, int *e2P, int n1, int n2, int e1
            unlockNodes(gotlocks1, lockw, 0, lockw, size);
            free(gotlocks1);
            free(lockw);
-           //CkPrintf("Possibly a livelock in edge_contract_help\n");
+           CkPrintf("Possibly a livelock in edge_contract_help\n");
            delete [] eConn;
            if(nesize!=0) delete[] nbrElems;
            free(gotlocks);
index b458914a6785bdd1d882e16cb0e18597d58d8ac8..f34f76887250d8ab38f0d8c31de92d1c49157afd 100644 (file)
@@ -1076,7 +1076,7 @@ void FEM_remove_element_local(FEM_Mesh *m, int element, int etype){
     The idxl info will be deleted in Purge, that way shared 
     chunks could still refer to it, needed for solution transfer
 */
-int FEM_remove_element(FEM_Mesh *m, int elementid, int elemtype, int permanent){
+int FEM_remove_element(FEM_Mesh *m, int elementid, int elemtype, int permanent, bool aggressive_node_removal){
   //CkAssert(elementid != -1);
 #ifdef DEBUG_2
   CkPrintf("removeElement, line %d\n", __LINE__);
@@ -1097,6 +1097,7 @@ int FEM_remove_element(FEM_Mesh *m, int elementid, int elemtype, int permanent){
     rm->elementid = sharedIdx;
     rm->elemtype = elemtype;
     rm->permanent = permanent;
+    rm->aggressive_node_removal = aggressive_node_removal;
     meshMod[remoteChunk].removeElementRemote(rm);
     // remove local ghost element, now done in purge
 #ifdef DEBUG_2
@@ -1470,23 +1471,26 @@ int FEM_remove_element(FEM_Mesh *m, int elementid, int elemtype, int permanent){
 #endif
 // This code removes nodes that we still need in cases involving acquisition/flip for
 // CPSD. I don't really follow the logic here, so I'm just cutting this bit out for now
-#ifndef CPSD
-                 if(!irec1) {
-                   if(!losinglocal) {
-                     FEM_remove_node_local(m,nodes[j]);
-                   }
-                   else {
-                     deleteNodeLater = nodes[j];
-                   }
+//#ifndef CPSD
+          CkPrintf("removeElement, aggressive_node_removal: %d\n", aggressive_node_removal);
+          if (aggressive_node_removal) {
+              if(!irec1) {
+                  if(!losinglocal) {
+                      FEM_remove_node_local(m,nodes[j]);
+                  }
+                  else {
+                      deleteNodeLater = nodes[j];
+                  }
 #ifdef DEBUG_2
-  CkPrintf("removeElement, line %d\n", __LINE__);
-#endif
-                   //if losing a local node, then can remove it only after its attributes 
-                   //have been copied, since this is the only copy.. 
-                   //(as no other chunk has this node as local/shared)
-                   //so we'll need to delete the ghostsend idxl entry and the node later
-                 }
+                  CkPrintf("removeElement, line %d\n", __LINE__);
 #endif
+                  //if losing a local node, then can remove it only after its attributes 
+                  //have been copied, since this is the only copy.. 
+                  //(as no other chunk has this node as local/shared)
+                  //so we'll need to delete the ghostsend idxl entry and the node later
+              }
+          }
+//#endif
                }
              }
              if(numElems!=0) delete[] elems;
@@ -2415,7 +2419,7 @@ void femMeshModify::removeGhostElem(removeGhostElemMsg *fm) {
 
 void femMeshModify::removeElementRemote(removeElemMsg *fm) {
   CtvAccess(_curTCharm) = tc;
-  fmUtil->removeElemRemote(fmMesh, fm->chk, fm->elementid, fm->elemtype, fm->permanent);
+  fmUtil->removeElemRemote(fmMesh, fm->chk, fm->elementid, fm->elemtype, fm->permanent, fm->aggressive_node_removal);
   delete fm;
   return;
 }
index bf41f1629314c7dbd1873951d25795d6a0e7eb47..8ebbbc626e9c8d4e07ea0ff2e988e8d41db61194 100644 (file)
@@ -626,10 +626,11 @@ void FEM_MUtil::addGhostElementRemote(FEM_Mesh *m, int chk, int elemType, int *i
 
 /** Remove this element, permanent specifies if this element is permanently removed
  */
-void FEM_MUtil::removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent) {
+void FEM_MUtil::removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent,
+        bool aggressive_node_removal) {
   const IDXL_List ll = m->elem[elemtype].ghostSend.addList(chk);
   int localIdx = ll[elementid];
-  FEM_remove_element(m, localIdx, elemtype, permanent);
+  FEM_remove_element(m, localIdx, elemtype, permanent, aggressive_node_removal);
   return;
 }
 
@@ -913,13 +914,14 @@ void FEM_MUtil::addToSharedList(FEM_Mesh *m, int fromChk, int sharedIdx) {
     Add a new element with this connectivity ON THIS CHUNK
     Update the lock on any node which changes owners
  */
-int FEM_MUtil::eatIntoElement(int localIdx) {
+int FEM_MUtil::eatIntoElement(int localIdx, bool aggressive_node_removal) {
   CkAssert(FEM_Is_ghost_index(localIdx));
   int nodesPerEl = mmod->fmMesh->elem[0].getConn().width();
   int *adjnodes = new int[nodesPerEl];
   int *oldnodes = new int[nodesPerEl];
   int elemtype = 0;
   mmod->fmMesh->e2n_getAll(localIdx, adjnodes, elemtype);
+  CkPrintf("Chunk %d eating elem %d(%d,%d,%d) %d\n",idx,localIdx,adjnodes[0],adjnodes[1],adjnodes[2], aggressive_node_removal);
 #ifdef DEBUG_1
   CkPrintf("Chunk %d eating elem %d(%d,%d,%d)\n",idx,localIdx,adjnodes[0],adjnodes[1],adjnodes[2]);
 #endif
@@ -937,7 +939,7 @@ int FEM_MUtil::eatIntoElement(int localIdx) {
 #ifdef DEBUG_1
   CkPrintf("eatIntoElement::remove\n");
 #endif
-  FEM_remove_element(mmod->fmMesh,localIdx,elemtype,idx);
+  FEM_remove_element(mmod->fmMesh,localIdx,elemtype,idx,aggressive_node_removal);
 #ifdef DEBUG_1
   CkPrintf("eatIntoElement::done removing\n");
 #endif
@@ -965,14 +967,14 @@ int FEM_MUtil::eatIntoElement(int localIdx) {
 #endif
   copyElemData(0,localIdx,newEl); //special copy across chunk
   FEM_purge_element(mmod->fmMesh,localIdx,elemtype);
-#ifndef CPSD
+//#ifndef CPSD
   for(int i=0; i<nodesPerEl; i++) {
     if(adjnodes[i]!=oldnodes[i]) {
       //correct the lock
       FEM_Modify_LockUpdate(mmod->fmMesh,adjnodes[i]);
     }
   }
-#endif
+//#endif
   free(adjnodes);
   free(oldnodes);
   return newEl;
@@ -1380,302 +1382,302 @@ void FEM_MUtil::FEM_Print_coords(FEM_Mesh *m, int nodeid) {
     Validity of nodes/elements
 */
 void FEM_MUtil::StructureTest(FEM_Mesh *m) {
-  int noNodes = m->node.size();
-  int noEle = m->elem[0].size();
-  int noGhostEle = m->elem[0].ghost->size();
-  int noGhostNodes = m->node.ghost->size();
-  int wdt = m->elem[0].getConn().width();
-  int *e2n = (int*)malloc(wdt*sizeof(int));
-  for(int i=0; i<noEle; i++) {
-    if(m->elem[0].is_valid(i)) {
-      m->e2n_getAll(i,e2n,0);
-      //must have all different connections
-      if(e2n[0]==e2n[1] || e2n[1]==e2n[2] || e2n[2]==e2n[0]) {
-       CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",i,e2n[0],e2n[1],e2n[2]);
-       CkAssert(false);
-      }
-      //local elem must have all local node connectivity 
-      if(e2n[0]<0 || e2n[1]<0 || e2n[2]<0) {
-       CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",i,e2n[0],e2n[1],e2n[2]);
-       CkAssert(false);
-      }
-      for(int j=0; j<3; j++) {
-       //all nodes must be valid and all local elems should have local connectivity
-       CkAssert(m->node.is_valid(e2n[j]));
-      }
-      int e2e[3];
-      m->e2e_getAll(i,e2e,0);
-      if((e2e[0]==e2e[1]) && (e2e[0]==e2e[2]) && (e2e[0]!=-1)) {
-       CkPrintf("ERROR: element %d, has e2e (%d,%d,%d)\n",i,e2e[0],e2e[1],e2e[2]);
-       CkAssert(false);
-      }
-    }
-  }
-  for(int i=0; i<noGhostEle; i++) {
-    if(m->elem[0].ghost->is_valid(i)) {
-      int ghostIndex = FEM_To_ghost_index(i);
-      m->e2n_getAll(ghostIndex,e2n,0);
-      if(e2n[0]==e2n[1] || e2n[1]==e2n[2] || e2n[2]==e2n[0]) {
-       //must have all different connections
-       CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",ghostIndex,e2n[0],e2n[1],e2n[2]);
-       CkAssert(false);
-      }
-      if(!(e2n[0]>=0 || e2n[1]>=0 || e2n[2]>=0)) {
-       //must have at least one local node
-       CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",ghostIndex,e2n[0],e2n[1],e2n[2]);
-       CkAssert(false);
-      }
-      for(int j=0; j<3; j++) {
-       //all nodes must be valid
-       if(e2n[j]>=0) CkAssert(m->node.is_valid(e2n[j]));
-       else CkAssert(m->node.ghost->is_valid(FEM_From_ghost_index(e2n[j])));
-      }
-      int e2e[3];
-      m->e2e_getAll(ghostIndex,e2e,0);
-      if((e2e[0]==e2e[1]) && (e2e[0]==e2e[2]) && (e2e[0]!=-1)) {
-       CkPrintf("ERROR: element %d, has e2e (%d,%d,%d)\n",i,e2e[0],e2e[1],e2e[2]);
-       CkAssert(false);
-      }
+    int noNodes = m->node.size();
+    int noEle = m->elem[0].size();
+    int noGhostEle = m->elem[0].ghost->size();
+    int noGhostNodes = m->node.ghost->size();
+    int wdt = m->elem[0].getConn().width();
+    int *e2n = (int*)malloc(wdt*sizeof(int));
+    for(int i=0; i<noEle; i++) {
+        if(m->elem[0].is_valid(i)) {
+            m->e2n_getAll(i,e2n,0);
+            //must have all different connections
+            if(e2n[0]==e2n[1] || e2n[1]==e2n[2] || e2n[2]==e2n[0]) {
+                CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",i,e2n[0],e2n[1],e2n[2]);
+                CkAssert(false);
+            }
+            //local elem must have all local node connectivity 
+            if(e2n[0]<0 || e2n[1]<0 || e2n[2]<0) {
+                CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",i,e2n[0],e2n[1],e2n[2]);
+                CkAssert(false);
+            }
+            for(int j=0; j<3; j++) {
+                //all nodes must be valid and all local elems should have local connectivity
+                CkAssert(m->node.is_valid(e2n[j]));
+            }
+            int e2e[3];
+            m->e2e_getAll(i,e2e,0);
+            if((e2e[0]==e2e[1]) && (e2e[0]==e2e[2]) && (e2e[0]!=-1)) {
+                CkPrintf("ERROR: element %d, has e2e (%d,%d,%d)\n",i,e2e[0],e2e[1],e2e[2]);
+                CkAssert(false);
+            }
+        }
     }
-  }
-  int *n2e, n2esize=0;
-  int *n2n, n2nsize=0;
-  for(int i=0; i<noNodes; i++) {
-    if(m->node.is_valid(i)) {
-      m->n2e_getAll(i,n2e,n2esize);
-      m->n2n_getAll(i,n2n,n2nsize);
-      if(n2esize>n2nsize) {
-       FEM_Print_coords(m,i);
-       FEM_Print_n2e(m,i);
-       FEM_Print_n2n(m,i);
-       CkPrintf("ERROR: local node %d, with inconsistent adjacency list\n",i);
-       CkAssert(false);
-      }
-    }
-    else {
-      continue;
-    }
-    if(n2esize > 0) {
-      for(int j=0; j<n2esize; j++) {
-       CkAssert(n2e[j]!=-1);
-       if(FEM_Is_ghost_index(n2e[j])) CkAssert(m->elem[0].ghost->is_valid(FEM_From_ghost_index(n2e[j]))==1);
-       else CkAssert(m->elem[0].is_valid(n2e[j])==1);
-      }
-      m->e2n_getAll(n2e[0],e2n,0);
-      //any local/shared node should have at least one local element connected to it
-      bool done = false;
-      for(int j=0; j<n2esize; j++) {
-       if(n2e[j] >= 0) {
-         done = true; 
-         break;
-       }
-      }
-      if(!done) {
-       FEM_Print_coords(m,i);
-       FEM_Print_n2e(m,i);
-       FEM_Print_n2n(m,i);
-       CkPrintf("ERROR: isolated local node %d, with no local element connectivity\n",i);
-       CkAssert(false);
-      }
-      //ensure that there is a cloud of connectivity, no disconnected elements, other than boundaries
-      int testnode = i;
-      int startnode = (e2n[0]==testnode) ? e2n[1] : e2n[0];
-      int othernode = (e2n[2]==testnode) ? e2n[1] : e2n[2];
-      int previousnode = startnode;
-      int nextnode = -1;
-      int numdeadends = 0;
-      int numunused = n2esize-1;
-      n2e[0] = -1;
-      for(int j=0; j<n2esize-1; j++) {
-       nextnode = -1;
-       for(int k=1; k<n2esize; k++) {
-         if(n2e[k]==-1) continue;
-         m->e2n_getAll(n2e[k],e2n,0);
-         if(e2n[0]==previousnode || e2n[1]==previousnode || e2n[2]==previousnode) {
-           nextnode = (e2n[0]==previousnode) ? ((e2n[1]==testnode)? e2n[2]:e2n[1]) : ((e2n[1]==previousnode)? ((e2n[0]==testnode)? e2n[2]:e2n[0]):((e2n[1]==testnode)? e2n[0]:e2n[1]));
-           previousnode = nextnode;
-           n2e[k] = -1;
-           numunused--;
-         }
-       }
-       if(nextnode==othernode && othernode!=-1) {
-         //it has reached a full circle
-         break;
-       }
-       else if(nextnode==-1) {
-         //this is one edge, start travelling along the other end
-         numdeadends++;
-         previousnode = othernode;
-         othernode = -1;
-       }
-       if(numdeadends>2 && numunused!=0) {
-         FEM_Print_coords(m,i);
-         FEM_Print_n2e(m,i);
-         FEM_Print_n2n(m,i);
-         CkPrintf("ERROR: cloud connectivity of node %d is discontinuous\n",i);
-         CkAssert(false);
-       }
-      }
-      if(n2esize>0) delete[] n2e; n2esize=0;
-      //reconstruct n2n from n2e & e2n
-      int n2n1size = 0;
-      int *n2n1 = (int*)malloc(n2nsize*sizeof(int));
-      int n2n1Count = 0;
-      m->n2e_getAll(i,n2e,n2esize);
-      for(int j=0; j<n2esize; j++) {
-       CkAssert(n2e[j]!=-1);
-       //each of these elems should have me in its e2n
-       int e2n1[3];
-       m->e2n_getAll(n2e[j],e2n1,0);
-       if(e2n1[0]!=i && e2n1[1]!=i && e2n1[2]!=i) {
-         FEM_Print_coords(m,i);
-         FEM_Print_n2e(m,i);
-         FEM_Print_e2n(m,n2e[j]);
-         CkPrintf("ERROR: ghost elem %d & ghost node %d have inconsistent adjacency list\n",n2e[j],i);
-         CkAssert(false);
-       }
-       for(int k=0; k<3;k++) {
-         if(e2n1[k] == i) continue;
-         bool flag1 = true;
-         for(int l=0; l<n2n1Count; l++) {
-           if(e2n1[k] == n2n1[l]) flag1 = false;
-         }
-         if(flag1 && n2n1Count<n2nsize) { //this is not in the list
-           n2n1[n2n1Count] = e2n1[k];
-           n2n1Count++;
-         }
-       }
-      }
-      //verify if n2n1 has the same nodes as n2n
-      bool flag2 = true;
-      if(n2n1Count!=n2nsize) flag2 = false;
-      for(int j=0; j<n2n1Count; j++) {
-       bool flag1 = false;
-       for(int k=0; k<n2nsize; k++) {
-         if(n2n[k]==n2n1[j]) flag1 = true;
-       }
-       if(!flag1) {
-         flag2 = false;
-         break;
-       }
-      }
-      if(!flag2) {
-       FEM_Print_coords(m,i);
-       FEM_Print_n2n(m,i);
-       for(int l=0; l<n2esize; l++) FEM_Print_e2n(m,n2e[l]);
-       CkPrintf("ERROR: ghost node %d has inconsistent adjacency list\n",i);
-       CkAssert(false);
-      }
-      
-      delete [] n2n1;
-      if(n2esize>0) delete [] n2e; n2esize=0;
-      if(n2nsize>0) delete [] n2n; n2nsize=0;
+    for(int i=0; i<noGhostEle; i++) {
+        if(m->elem[0].ghost->is_valid(i)) {
+            int ghostIndex = FEM_To_ghost_index(i);
+            m->e2n_getAll(ghostIndex,e2n,0);
+            if(e2n[0]==e2n[1] || e2n[1]==e2n[2] || e2n[2]==e2n[0]) {
+                //must have all different connections
+                CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",ghostIndex,e2n[0],e2n[1],e2n[2]);
+                CkAssert(false);
+            }
+            if(!(e2n[0]>=0 || e2n[1]>=0 || e2n[2]>=0)) {
+                //must have at least one local node
+                CkPrintf("ERROR: element %d, has connectivity (%d,%d,%d)\n",ghostIndex,e2n[0],e2n[1],e2n[2]);
+                CkAssert(false);
+            }
+            for(int j=0; j<3; j++) {
+                //all nodes must be valid
+                if(e2n[j]>=0) CkAssert(m->node.is_valid(e2n[j]));
+                else CkAssert(m->node.ghost->is_valid(FEM_From_ghost_index(e2n[j])));
+            }
+            int e2e[3];
+            m->e2e_getAll(ghostIndex,e2e,0);
+            if((e2e[0]==e2e[1]) && (e2e[0]==e2e[2]) && (e2e[0]!=-1)) {
+                CkPrintf("ERROR: element %d, has e2e (%d,%d,%d)\n",i,e2e[0],e2e[1],e2e[2]);
+                CkAssert(false);
+            }
+        }
     }
-  }
-  if(n2esize>0) delete [] n2e; n2esize=0;
-  if(n2nsize>0) delete [] n2n; n2nsize=0;
-  for(int i=0; i<noGhostNodes; i++) {
-    int ghostidx = FEM_To_ghost_index(i);
-    if(m->node.ghost->is_valid(i)) {
-      m->n2e_getAll(ghostidx,n2e,n2esize);
-      m->n2n_getAll(ghostidx,n2n,n2nsize);
-      bool done = false;
-      for(int k=0;k<n2nsize;k++) {
-       if(n2n[k]>=0) {
-         done = true;
-         break;
-       }
-      }
-      if(n2esize>n2nsize || !done) {
-       FEM_Print_coords(m,ghostidx);
-       FEM_Print_n2e(m,ghostidx);
-       FEM_Print_n2n(m,ghostidx);
-       CkPrintf("ERROR: ghost node %d, with inconsistent adjacency list\n",ghostidx);
-       CkAssert(false);
-      }
-      if(n2esize > 0) {
-       //reconstruct n2n from n2e & e2n
-       int n2n1size = 0;
-       int *n2n1 = (int*)malloc(n2nsize*sizeof(int));
-       int n2n1Count = 0;
-       for(int j=0; j<n2esize; j++) {
-         CkAssert(n2e[j]!=-1);
-         if(FEM_Is_ghost_index(n2e[j])) {
-           CkAssert(m->elem[0].ghost->is_valid(FEM_From_ghost_index(n2e[j]))==1);
-         }
-         else {
-           CkAssert(m->elem[0].is_valid(n2e[j])==1);
-         }
-         //each of these elems should have me in its e2n
-         int e2n1[3];
-         m->e2n_getAll(n2e[j],e2n1,0);
-         if(e2n1[0]!=ghostidx && e2n1[1]!=ghostidx && e2n1[2]!=ghostidx) {
-           FEM_Print_coords(m,ghostidx);
-           FEM_Print_n2e(m,ghostidx);
-           FEM_Print_e2n(m,n2e[j]);
-           CkPrintf("ERROR: ghost elem %d & ghost node %d have inconsistent adjacency list\n",n2e[j],ghostidx);
-           CkAssert(false);
-         }
-         for(int k=0; k<3;k++) {
-           if(e2n1[k] == ghostidx) continue;
-           bool flag1 = true;
-           for(int l=0; l<n2n1Count; l++) {
-             if(e2n1[k] == n2n1[l]) flag1 = false;
-           }
-           if(flag1 && n2n1Count<n2nsize) { //this is not in the list
-             n2n1[n2n1Count] = e2n1[k];
-             n2n1Count++;
-           }
-         }
-       }
-       //verify if n2n1 has the same nodes as n2n
-       bool flag2 = true;
-       if(n2n1Count!=n2nsize) flag2 = false;
-       for(int j=0; j<n2n1Count; j++) {
-         bool flag1 = false;
-         for(int k=0; k<n2nsize; k++) {
-           if(n2n[k]==n2n1[j]) flag1 = true;
-         }
-         if(!flag1) {
-           flag2 = false;
-           break;
-         }
-       }
-       if(!flag2) {
-         FEM_Print_coords(m,ghostidx);
-         FEM_Print_n2n(m,ghostidx);
-         for(int l=0; l<n2esize; l++) FEM_Print_e2n(m,n2e[l]);
-         CkPrintf("ERROR: ghost node %d has inconsistent adjacency list\n",ghostidx);
-         CkAssert(false);
-       }
-       delete[] n2n1;
-       delete[] n2e;
-      }
-      if(n2nsize > 0) {
-       for(int j=0; j<n2nsize; j++) {
-         CkAssert(n2n[j]!=-1);
-       }
-       delete[] n2n;
-      }
-      //verify that it is coming correctly from all chunks as a ghost
-      const IDXL_Rec *irec = m->node.ghost->ghostRecv.getRec(i);
-      //go to a chunk which owns it & verify that these are the only chunks that own it
-      int remChk = irec->getChk(0);
-      int sharedIdx = irec->getIdx(0);
-      int numsh = irec->getShared();
-      verifyghostsendMsg *vmsg = new(numsh,0)verifyghostsendMsg();
-      vmsg->fromChk = idx;
-      vmsg->sharedIdx = sharedIdx;
-      vmsg->numchks = numsh;
-      for(int k=0; k<numsh; k++) vmsg->chunks[k] = irec->getChk(k);
-      meshMod[remChk].verifyghostsend(vmsg);
+    int *n2e, n2esize=0;
+    int *n2n, n2nsize=0;
+    for(int i=0; i<noNodes; i++) {
+        if(m->node.is_valid(i)) {
+            m->n2e_getAll(i,n2e,n2esize);
+            m->n2n_getAll(i,n2n,n2nsize);
+            if(n2esize>n2nsize) {
+                FEM_Print_coords(m,i);
+                FEM_Print_n2e(m,i);
+                FEM_Print_n2n(m,i);
+                CkPrintf("ERROR: local node %d, with inconsistent adjacency list\n",i);
+                CkAssert(false);
+            }
+        }
+        else {
+            continue;
+        }
+        if(n2esize > 0) {
+            for(int j=0; j<n2esize; j++) {
+                CkAssert(n2e[j]!=-1);
+                if(FEM_Is_ghost_index(n2e[j])) CkAssert(m->elem[0].ghost->is_valid(FEM_From_ghost_index(n2e[j]))==1);
+                else CkAssert(m->elem[0].is_valid(n2e[j])==1);
+            }
+            m->e2n_getAll(n2e[0],e2n,0);
+            //any local/shared node should have at least one local element connected to it
+            bool done = false;
+            for(int j=0; j<n2esize; j++) {
+                if(n2e[j] >= 0) {
+                    done = true; 
+                    break;
+                }
+            }
+            if(!done) {
+                FEM_Print_coords(m,i);
+                FEM_Print_n2e(m,i);
+                FEM_Print_n2n(m,i);
+                CkPrintf("ERROR: isolated local node %d, with no local element connectivity\n",i);
+                CkAssert(false);
+            }
+            //ensure that there is a cloud of connectivity, no disconnected elements, other than boundaries
+            int testnode = i;
+            int startnode = (e2n[0]==testnode) ? e2n[1] : e2n[0];
+            int othernode = (e2n[2]==testnode) ? e2n[1] : e2n[2];
+            int previousnode = startnode;
+            int nextnode = -1;
+            int numdeadends = 0;
+            int numunused = n2esize-1;
+            n2e[0] = -1;
+            for(int j=0; j<n2esize-1; j++) {
+                nextnode = -1;
+                for(int k=1; k<n2esize; k++) {
+                    if(n2e[k]==-1) continue;
+                    m->e2n_getAll(n2e[k],e2n,0);
+                    if(e2n[0]==previousnode || e2n[1]==previousnode || e2n[2]==previousnode) {
+                        nextnode = (e2n[0]==previousnode) ? ((e2n[1]==testnode)? e2n[2]:e2n[1]) : ((e2n[1]==previousnode)? ((e2n[0]==testnode)? e2n[2]:e2n[0]):((e2n[1]==testnode)? e2n[0]:e2n[1]));
+                        previousnode = nextnode;
+                        n2e[k] = -1;
+                        numunused--;
+                    }
+                }
+                if(nextnode==othernode && othernode!=-1) {
+                    //it has reached a full circle
+                    break;
+                }
+                else if(nextnode==-1) {
+                    //this is one edge, start travelling along the other end
+                    numdeadends++;
+                    previousnode = othernode;
+                    othernode = -1;
+                }
+                if(numdeadends>2 && numunused!=0) {
+                    FEM_Print_coords(m,i);
+                    FEM_Print_n2e(m,i);
+                    FEM_Print_n2n(m,i);
+                    CkPrintf("ERROR: cloud connectivity of node %d is discontinuous\n",i);
+                    CkAssert(false);
+                }
+            }
+            if(n2esize>0) delete[] n2e; n2esize=0;
+            //reconstruct n2n from n2e & e2n
+            int n2n1size = 0;
+            int *n2n1 = (int*)malloc(n2nsize*sizeof(int));
+            int n2n1Count = 0;
+            m->n2e_getAll(i,n2e,n2esize);
+            for(int j=0; j<n2esize; j++) {
+                CkAssert(n2e[j]!=-1);
+                //each of these elems should have me in its e2n
+                int e2n1[3];
+                m->e2n_getAll(n2e[j],e2n1,0);
+                if(e2n1[0]!=i && e2n1[1]!=i && e2n1[2]!=i) {
+                    FEM_Print_coords(m,i);
+                    FEM_Print_n2e(m,i);
+                    FEM_Print_e2n(m,n2e[j]);
+                    CkPrintf("ERROR: ghost elem %d & ghost node %d have inconsistent adjacency list\n",n2e[j],i);
+                    CkAssert(false);
+                }
+                for(int k=0; k<3;k++) {
+                    if(e2n1[k] == i) continue;
+                    bool flag1 = true;
+                    for(int l=0; l<n2n1Count; l++) {
+                        if(e2n1[k] == n2n1[l]) flag1 = false;
+                    }
+                    if(flag1 && n2n1Count<n2nsize) { //this is not in the list
+                        n2n1[n2n1Count] = e2n1[k];
+                        n2n1Count++;
+                    }
+                }
+            }
+            //verify if n2n1 has the same nodes as n2n
+            bool flag2 = true;
+            if(n2n1Count!=n2nsize) flag2 = false;
+            for(int j=0; j<n2n1Count; j++) {
+                bool flag1 = false;
+                for(int k=0; k<n2nsize; k++) {
+                    if(n2n[k]==n2n1[j]) flag1 = true;
+                }
+                if(!flag1) {
+                    flag2 = false;
+                    break;
+                }
+            }
+            if(!flag2) {
+                FEM_Print_coords(m,i);
+                FEM_Print_n2n(m,i);
+                for(int l=0; l<n2esize; l++) FEM_Print_e2n(m,n2e[l]);
+                CkPrintf("ERROR: ghost node %d has inconsistent adjacency list\n",i);
+                CkAssert(false);
+            }
+
+            delete [] n2n1;
+            if(n2esize>0) delete [] n2e; n2esize=0;
+            if(n2nsize>0) delete [] n2n; n2nsize=0;
+        }
     }
-    else {
-      continue;
+    if(n2esize>0) delete [] n2e; n2esize=0;
+    if(n2nsize>0) delete [] n2n; n2nsize=0;
+    for(int i=0; i<noGhostNodes; i++) {
+        int ghostidx = FEM_To_ghost_index(i);
+        if(m->node.ghost->is_valid(i)) {
+            m->n2e_getAll(ghostidx,n2e,n2esize);
+            m->n2n_getAll(ghostidx,n2n,n2nsize);
+            bool done = false;
+            for(int k=0;k<n2nsize;k++) {
+                if(n2n[k]>=0) {
+                    done = true;
+                    break;
+                }
+            }
+            if(n2esize>n2nsize || !done) {
+                FEM_Print_coords(m,ghostidx);
+                FEM_Print_n2e(m,ghostidx);
+                FEM_Print_n2n(m,ghostidx);
+                CkPrintf("ERROR: ghost node %d, with inconsistent adjacency list\n",ghostidx);
+                CkAssert(false);
+            }
+            if(n2esize > 0) {
+                //reconstruct n2n from n2e & e2n
+                int n2n1size = 0;
+                int *n2n1 = (int*)malloc(n2nsize*sizeof(int));
+                int n2n1Count = 0;
+                for(int j=0; j<n2esize; j++) {
+                    CkAssert(n2e[j]!=-1);
+                    if(FEM_Is_ghost_index(n2e[j])) {
+                        CkAssert(m->elem[0].ghost->is_valid(FEM_From_ghost_index(n2e[j]))==1);
+                    }
+                    else {
+                        CkAssert(m->elem[0].is_valid(n2e[j])==1);
+                    }
+                    //each of these elems should have me in its e2n
+                    int e2n1[3];
+                    m->e2n_getAll(n2e[j],e2n1,0);
+                    if(e2n1[0]!=ghostidx && e2n1[1]!=ghostidx && e2n1[2]!=ghostidx) {
+                        FEM_Print_coords(m,ghostidx);
+                        FEM_Print_n2e(m,ghostidx);
+                        FEM_Print_e2n(m,n2e[j]);
+                        CkPrintf("ERROR: ghost elem %d & ghost node %d have inconsistent adjacency list\n",n2e[j],ghostidx);
+                        CkAssert(false);
+                    }
+                    for(int k=0; k<3;k++) {
+                        if(e2n1[k] == ghostidx) continue;
+                        bool flag1 = true;
+                        for(int l=0; l<n2n1Count; l++) {
+                            if(e2n1[k] == n2n1[l]) flag1 = false;
+                        }
+                        if(flag1 && n2n1Count<n2nsize) { //this is not in the list
+                            n2n1[n2n1Count] = e2n1[k];
+                            n2n1Count++;
+                        }
+                    }
+                }
+                //verify if n2n1 has the same nodes as n2n
+                bool flag2 = true;
+                if(n2n1Count!=n2nsize) flag2 = false;
+                for(int j=0; j<n2n1Count; j++) {
+                    bool flag1 = false;
+                    for(int k=0; k<n2nsize; k++) {
+                        if(n2n[k]==n2n1[j]) flag1 = true;
+                    }
+                    if(!flag1) {
+                        flag2 = false;
+                        break;
+                    }
+                }
+                if(!flag2) {
+                    FEM_Print_coords(m,ghostidx);
+                    FEM_Print_n2n(m,ghostidx);
+                    for(int l=0; l<n2esize; l++) FEM_Print_e2n(m,n2e[l]);
+                    CkPrintf("ERROR: ghost node %d has inconsistent adjacency list\n",ghostidx);
+                    CkAssert(false);
+                }
+                delete[] n2n1;
+                delete[] n2e;
+            }
+            if(n2nsize > 0) {
+                for(int j=0; j<n2nsize; j++) {
+                    CkAssert(n2n[j]!=-1);
+                }
+                delete[] n2n;
+            }
+            //verify that it is coming correctly from all chunks as a ghost
+            const IDXL_Rec *irec = m->node.ghost->ghostRecv.getRec(i);
+            //go to a chunk which owns it & verify that these are the only chunks that own it
+            int remChk = irec->getChk(0);
+            int sharedIdx = irec->getIdx(0);
+            int numsh = irec->getShared();
+            verifyghostsendMsg *vmsg = new(numsh,0)verifyghostsendMsg();
+            vmsg->fromChk = idx;
+            vmsg->sharedIdx = sharedIdx;
+            vmsg->numchks = numsh;
+            for(int k=0; k<numsh; k++) vmsg->chunks[k] = irec->getChk(k);
+            meshMod[remChk].verifyghostsend(vmsg);
+        }
+        else {
+            continue;
+        }
     }
-  }
-  free(e2n);
-  return;
+    free(e2n);
+    return;
 }
 
 /** The area test verifies that the area of no element is less than the SLIVERAREA
@@ -1792,7 +1794,10 @@ int FEM_MUtil::residualLockTest(FEM_Mesh *m) {
   int noNodes = m->node.size();
   for(int i=0; i<noNodes; i++) {
     if(m->node.is_valid(i)) {
-      CkAssert(!mmod->fmLockN[i].haslocks());
+      if (mmod->fmLockN[i].haslocks()) {
+          CkPrintf("[%d] Node %d has a residual lock\n", FEM_My_partition(), i);
+          CkAssert(false);
+      }
     }
   }
   for(int i=0; i<mmod->numChunks; i++) {
@@ -1800,3 +1805,22 @@ int FEM_MUtil::residualLockTest(FEM_Mesh *m) {
   }
   return 1;
 }
+
+/**
+ * Remove all extant node locks in the mesh. Probably only useful for debugging
+ * and very special cases.
+ */
+void FEM_MUtil::unlockAll(FEM_Mesh *m) {
+  int noNodes = m->node.size();
+  for(int i=0; i<noNodes; i++) {
+    if(m->node.is_valid(i)) {
+      if (mmod->fmLockN[i].haslocks()) {
+          mmod->fmLockN[i].reset(i, mmod);
+      }
+    }
+  }
+  for(int i=0; i<mmod->numChunks; i++) {
+    mmod->fmIdxlLock[i] = false;
+  }
+}
+