[MOAB-dev] r5579 - in MOAB/trunk/src: . io moab

tautges at mcs.anl.gov tautges at mcs.anl.gov
Thu Jun 14 13:25:14 CDT 2012


Author: tautges
Date: 2012-06-14 13:25:14 -0500 (Thu, 14 Jun 2012)
New Revision: 5579

Modified:
   MOAB/trunk/src/ScdInterface.cpp
   MOAB/trunk/src/io/ReadNC.cpp
   MOAB/trunk/src/moab/ScdInterface.hpp
Log:
Fix the scdpart test, and do periodic structured mesh correctly in parallel.
In serial, a periodic structured mesh is represented that way, i.e. in the i direction, 
the last column of i vertices corresponds to the first column.  In parallel, the local
mesh is non-periodic, but globally, for periodic in i, the last (global) column of i vertices on
the proc holding that column should be identified with the first (global) column of i vertices
on the proc holding that column.  This wasn't really the case before, so parallel structured
periodic mesh didn't work right.

For now, the fix is in the NC reader.  One could argue that this should really be implemented
in ScdInterface.  This commit also doesn't include a detailed test case that catches the
problem (though scdpart was failing before, but succeeds now).  So there may be more coming.


Index: src/io/ReadNC.cpp
===================================================================
--- src/io/ReadNC.cpp	(revision 5573)
+++ src/io/ReadNC.cpp	(working copy)
@@ -328,9 +328,14 @@
 {
   Range tmp_range;
   ScdBox *scd_box;
+  bool lperiodic_i = false, gperiodic_i = true;
+#ifdef USE_MPI
+    // if serial, use a locally-periodic representation, otherwise don't
+  if (myPcomm->proc_config().proc_size() == 1) lperiodic_i = true;
+#endif  
   ErrorCode rval = scdi->construct_box(HomCoord(lDims[0], lDims[1], (-1 != lDims[2] ? lDims[2] : 0), 1),
                                        HomCoord(lDims[3], lDims[4], (-1 != lDims[5] ? lDims[5] : 0), 1),
-                                       NULL, 0, scd_box, true);
+                                       NULL, 0, scd_box, lperiodic_i, false);
   ERRORR(rval, "Trouble creating scd vertex sequence.");
 
     // set the global box parameters
@@ -362,13 +367,14 @@
   rval = scd_box->get_coordinate_arrays(xc, yc, zc);
   ERRORR(rval, "Couldn't get vertex coordinate arrays.");
 
-  int i, j, k, il, jl, kl;
+  int i, j, k, il, jl, kl, itmp;
   int dil = lDims[3] - lDims[0] + 1;
   int djl = lDims[4] - lDims[1] + 1;
   int di = gDims[3] - gDims[0] + 1;
   int dj = gDims[4] - gDims[1] + 1;
   assert(dil == (int)ilVals.size() && djl == (int)jlVals.size() && 
          (-1 == lDims[2] || lDims[5]-lDims[2]+1 == (int)klVals.size()));
+#define INDEX(i,j,k) ()
   for (kl = lDims[2]; kl <= lDims[5]; kl++) {
     k = kl - lDims[2];
     for (jl = lDims[1]; jl <= lDims[4]; jl++) {
@@ -379,11 +385,13 @@
         xc[pos] = ilVals[i];
         yc[pos] = jlVals[j];
         zc[pos] = (-1 == lDims[2] ? 0.0 : klVals[k]);
-        *gid_data = (-1 != kl ? kl*di*dj : 0) + jl*di + il + 1;
+        itmp = (!lperiodic_i && gperiodic_i && il == gDims[3] ? gDims[0] : il);
+        *gid_data = (-1 != kl ? kl*di*dj : 0) + jl*di + itmp + 1;
         gid_data++;
       }
     }
   }
+#undef INDEX
 
 #ifndef NDEBUG
   int num_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) *
Index: src/ScdInterface.cpp
===================================================================
--- src/ScdInterface.cpp	(revision 5573)
+++ src/ScdInterface.cpp	(working copy)
@@ -759,6 +759,7 @@
 {
 #ifdef USE_MPI
   if (dijk[2] != 0) {
+      // for sqij, there is no k neighbor, ever
     pto = -1;
     return MB_SUCCESS;
   }
@@ -766,69 +767,82 @@
   pto = -1;
   bdy_ind[0] = bdy_ind[1] = -1;
   int pi(0), pj; // pi, pj: # procs in i, j directions
-              // guess pi
+              // get pi
   ErrorCode rval = compute_partition_sqij(np, nr, gdims, rdims, &pi);
   if (MB_SUCCESS != rval) return rval;
   pj = np / pi;
   assert(pi * pj == np);
   pto = -1;
-  if ((!(nr%pi) && -1 == dijk[0] && !periodic_i) ||  // left and not periodic
-      ((nr%pi) == pi-1 && 1 == dijk[0] && !periodic_i) ||  // right and not periodic
-      (!(nr/pi) && -1 == dijk[1] && !periodic_j) || // bottom and not periodic
-      (nr/pi == pj-1 && 1 == dijk[1] && !periodic_j))  // top and not periodic
+  bool top_i = 0, top_j = 0, bot_i = 0, bot_j = 0;
+  if (nr%pi == pi-1) top_i = 1;
+  if (nr/pi == pj-1) top_j = 1;
+  if (!(nr%pi)) bot_i = 1;
+  if (!(nr/pi)) bot_j = 1;
+  if ((!periodic_i && bot_i && -1 == dijk[0]) ||  // left and not periodic
+      (!periodic_i && top_i && 1 == dijk[0]) ||  // right and not periodic
+      (!periodic_j && bot_j && -1 == dijk[1]) || // bottom and not periodic
+      (!periodic_j && top_j && 1 == dijk[1]))  // top and not periodic
     return MB_SUCCESS;
   
   std::copy(ldims, ldims+6, facedims);
   pto = nr;
-  int dj = (gdims[4] - gdims[1]) / pj, jextra = (gdims[4] - gdims[1]) % dj,
-      di = (gdims[3] - gdims[0]) / pi, iextra = (gdims[3] - gdims[0]) % di;
+  int j = gdims[4] - gdims[1], dj = j / pj, jextra = (gdims[4] - gdims[1]) % dj,
+      i = gdims[3] - gdims[0], di = i / pi, iextra = (gdims[3] - gdims[0]) % di;
   
   if (0 != dijk[0]) {
-    pto = (pto + dijk[0]) % np;
+    pto = (pto + dijk[0] + np) % np;
     assert (pto >= 0 && pto < np);
     if (-1 == dijk[0]) {
       facedims[3] = facedims[0];
-      rdims[3] = rdims[0];
-      rdims[0] -= di;
+      rdims[3] = bot_i ? gdims[3] : rdims[0];
+      rdims[0] = rdims[3] - di;
       if (pto%pi < iextra) rdims[0]--;
     }
     else {
+      if (top_i) facedims[3] = gdims[0];
       facedims[0] = facedims[3];
-      rdims[0] = rdims[3];
-      rdims[3] += di;
+      rdims[0] = (top_i ? gdims[0] : rdims[3]);
+      rdims[3] = rdims[0] + di;
       if (pto%pi < iextra) rdims[3]++;
     }
   }
   if (0 != dijk[1]) {
-    pto = (pto + dijk[1]*pi) % np;
+    pto = (pto + dijk[1]*pi + np) % np;
     assert (pto >= 0 && pto < np);
     if (-1 == dijk[1]) {
       facedims[4] = facedims[1];
-      rdims[4] = rdims[1];
-      rdims[1] -= dj;
+      rdims[4] = (bot_j ? gdims[4] : rdims[1]);
+      rdims[1] = rdims[4] - dj;
       if (pto/pi < jextra) rdims[1]--;
     }
     else {
+      if (top_j) facedims[4] = gdims[1];
       facedims[1] = facedims[4];
-      rdims[1] = rdims[4];
-      rdims[4] += dj;
+      rdims[1] = (top_j ? gdims[1] : rdims[4]);
+      rdims[4] = rdims[1] + dj;
       if (pto/pi < jextra) rdims[4]++;
     }
   }
 
-  if (dijk[0] == -1 && periodic_i && ldims[0] == gdims[0]) bdy_ind[0] = 1;
-  if (dijk[1] == -1 && periodic_j && ldims[1] == gdims[1]) bdy_ind[1] = 1;
+  if (dijk[0] == -1 && periodic_i && bot_i) bdy_ind[0] = 1;
+  if (dijk[1] == -1 && periodic_j && bot_j) bdy_ind[1] = 1;
 
   assert(-1 == pto ||
          (rdims[0] >= gdims[0] && rdims[3] <= gdims[3] && 
           rdims[1] >= gdims[1] && rdims[4] <= gdims[4] && 
           rdims[2] >= gdims[2] && rdims[5] <= gdims[5] &&
-          facedims[0] >= rdims[0] && facedims[3] <= rdims[3] &&
-          facedims[1] >= rdims[1] && facedims[4] <= rdims[4] &&
-          facedims[2] >= rdims[2] && facedims[5] <= rdims[5] &&
-          facedims[0] >= ldims[0] && facedims[3] <= ldims[3] &&
-          facedims[1] >= ldims[1] && facedims[4] <= ldims[4] &&
-          facedims[2] >= ldims[2] && facedims[5] <= ldims[5]));
+          (facedims[0] >= rdims[0] || (periodic_i && rdims[3] == gdims[3] && facedims[0] == gdims[0])) && 
+          facedims[3] <= rdims[3] &&
+          (facedims[1] >= rdims[1]  || (periodic_j && rdims[4] == gdims[4] && facedims[1] == gdims[1])) && 
+          facedims[4] <= rdims[4] &&
+          facedims[2] >= rdims[2] && 
+          facedims[5] <= rdims[5] &&
+          facedims[0] >= ldims[0] && 
+          facedims[3] <= ldims[3] &&
+          facedims[1] >= ldims[1] && 
+          facedims[4] <= ldims[4] &&
+          facedims[2] >= ldims[2] && 
+          facedims[5] <= ldims[5]));
 
   return MB_SUCCESS;
 #else
@@ -1016,9 +1030,10 @@
                             i, j, k, pto, bdy_ind, ijkrem, ijkface);
         if (MB_SUCCESS != rval) return rval;
         if (-1 != pto) {
-          assert(std::find(procs.begin(), procs.end(), pto) == procs.end());
-          procs.push_back(pto);
-          offsets.push_back(shared_indices.size());
+          if (procs.empty() || pto != *procs.rbegin()) {
+            procs.push_back(pto);
+            offsets.push_back(shared_indices.size());
+          }
           rval = get_indices(bdy_ind, ldims, ijkrem, ijkface, shared_indices);
           if (MB_SUCCESS != rval) return rval;
         }
Index: src/moab/ScdInterface.hpp
===================================================================
--- src/moab/ScdInterface.hpp	(revision 5573)
+++ src/moab/ScdInterface.hpp	(working copy)
@@ -259,8 +259,8 @@
     /** Partitions the structured parametric space by seeking square ij partitions
      * For description of arguments, see ScdInterface::compute_partition.
      */
-  inline static ErrorCode compute_partition_sqij(int np, int nr, const int *gijk, int *lijk,
-                                                 int *pip = NULL);
+  inline static ErrorCode compute_partition_sqij(int np, int nr, const int *gijk, 
+                                                 int *lijk, int *pip = NULL);
   
     //! Compute a partition of structured parameter space
     /** Partitions the structured parametric space by seeking square jk partitions
@@ -766,8 +766,8 @@
   return MB_SUCCESS;
 }
 
-inline ErrorCode ScdInterface::compute_partition_sqij(int np, int nr, const int *gijk, int *lijk,
-                                                      int *pip) 
+inline ErrorCode ScdInterface::compute_partition_sqij(int np, int nr, const int *gijk, 
+                                                      int *lijk, int *pip) 
 {
     // square IxJ partition
 



Modified: MOAB/trunk/src/ScdInterface.cpp
===================================================================
--- MOAB/trunk/src/ScdInterface.cpp	2012-06-14 18:10:19 UTC (rev 5578)
+++ MOAB/trunk/src/ScdInterface.cpp	2012-06-14 18:25:14 UTC (rev 5579)
@@ -699,19 +699,20 @@
   std::copy(ldims, ldims+6, facedims);
   
   if (0 != dijk[1]) {
-    pto = (pto + dijk[1]*pk) % np;
+    pto = (pto + dijk[1]*pk + np) % np;
     assert (pto >= 0 && pto < np);
     int dj = (gdims[4] - gdims[1]) / pj, extra = (gdims[4] - gdims[1]) % pj;
     if (-1 == dijk[1]) {
       facedims[4] = facedims[1];
-      rdims[4] = rdims[1];
-      rdims[1] -= dj;
+      rdims[4] = (nr < pk ? gdims[4] : rdims[0]);
+      rdims[1] = rdims[4] - dj;
       if (pto < extra) rdims[1]--;
     }
     else {
+      if (nr > np-nk) facedims[4] = gdims[1];
       facedims[1] = facedims[4];
-      rdims[1] = rdims[4];
-      rdims[4] += dj;
+      rdims[1] = (nr > np-nk ? gdims[1] : rdims[4]);
+      rdims[4] = rdims[1] + dj;
       if (pto < extra) rdims[4]++;
     }
   }
@@ -738,12 +739,18 @@
          (rdims[0] >= gdims[0] && rdims[3] <= gdims[3] && 
           rdims[1] >= gdims[1] && rdims[4] <= gdims[4] && 
           rdims[2] >= gdims[2] && rdims[5] <= gdims[5] &&
-          facedims[0] >= rdims[0] && facedims[3] <= rdims[3] &&
-          facedims[1] >= rdims[1] && facedims[4] <= rdims[4] &&
-          facedims[2] >= rdims[2] && facedims[5] <= rdims[5] &&
-          facedims[0] >= ldims[0] && facedims[3] <= ldims[3] &&
-          facedims[1] >= ldims[1] && facedims[4] <= ldims[4] &&
-          facedims[2] >= ldims[2] && facedims[5] <= ldims[5]));
+          (facedims[0] >= rdims[0] || (periodic_i && rdims[3] == gdims[3] && facedims[0] == gdims[0])) && 
+          facedims[3] <= rdims[3] &&
+          (facedims[1] >= rdims[1]  || (periodic_j && rdims[4] == gdims[4] && facedims[1] == gdims[1])) && 
+          facedims[4] <= rdims[4] &&
+          facedims[2] >= rdims[2] && 
+          facedims[5] <= rdims[5] &&
+          facedims[0] >= ldims[0] && 
+          facedims[3] <= ldims[3] &&
+          facedims[1] >= ldims[1] && 
+          facedims[4] <= ldims[4] &&
+          facedims[2] >= ldims[2] && 


More information about the moab-dev mailing list