CarpetAdaptiveRegrid: Enforce proper nesting

Enforce proper nesting both for parents and children by looping over
all finer levels regridding them all at once. This completely destroys
efficiency...

darcs-hash:20050222174437-58c7f-d280bacafa2fb4c634613381b4c7134733348ec5.gz

1 files changed, 615 insertions, 468 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index 6923d9188..f99d44e84 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -105,11 +105,18 @@ namespace CarpetAdaptiveRegrid {
       if (cctkGH->cctk_iteration ==0) {
         return 0;
       }
-
-    }    
+      
+    }
 
     if (reflevel == maxreflevels - 1) return 0;
 
+    // Return if we want to regrid regularly, but not at this time
+    if (regrid_every > 0 && cctkGH->cctk_iteration != 0
+        && (cctkGH->cctk_iteration-1) % regrid_every != 0)
+    {
+      return 0;
+    }
+
 //     cout << "bbsss at start" << endl << bbsss << endl;
 //     cout << "obss at start" << endl << obss << endl;
 //     cout << "pss at start" << endl << pss << endl;
@@ -117,563 +124,703 @@ namespace CarpetAdaptiveRegrid {
     CCTK_INT do_recompose;
     do_recompose = 1;
 
-    // We assume that we are called in fine->coarse order
-    // So that on finer levels regridding has already been done
-
-    // So the full algorithm should look something like:
-
-    // Find how big the first bounding box should be on this level
-    //   Do this by finding min lower and max upper bounds of all bboxes
-    // Allocate box
-    // Fill errors from local arrays
-    // If grandchildren exist use their bboxes (expanded) to add to errors
-    // Reduce errors (MPI sum)
-    // Set errors to 0/1
-    // Define vectors of bboxes final (empty) and todo (contains bbox)
-    // Define vector of masks (contains error mask)
-    // Loop over all entries in todo:
-    //   Setup appropriate 1d array memory
-    //   Call fortran routine
-    //   If return is:
-    //     zero: add bbox to final
-    //     one:  add new bbox to todo and assoc mask to masklist
-    //     two:  add both new bboxs to todo and assoc masks to masklist
-    //   
-
-    vector<ibbox> bbs = local_bbsss.at(mglevel).at(reflevel);
+    CCTK_INT called_on_ml = mglevel;
+    CCTK_INT called_on_rl = reflevel;
+    CCTK_INT called_on_map = carpetGH.map;
     
-    stack<ibbox> final;
-
-    vector<vector<ibbox> > bbss = local_bbsss.at(0);
-      
-    //    ivect low = 100000, upp = -1;
-
-    bool did_regrid = false;
+    CCTK_INT finest_current_rl = local_bbsss.at(0).size();
+    finest_current_rl = min(finest_current_rl, maxreflevels - 1);
 
-    rvect physical_min, physical_max;
-    rvect interior_min, interior_max;
-    rvect exterior_min, exterior_max;
-    rvect base_spacing;
-    int ierr = GetDomainSpecification
-      (dim, &physical_min[0], &physical_max[0],
-       &interior_min[0], &interior_max[0],
-       &exterior_min[0], &exterior_max[0], &base_spacing[0]);
-    assert (!ierr);
-
-    for ( vector<ibbox>::const_iterator bbi = bbs.begin(); 
-          bbi != bbs.end();
-          ++bbi) 
-    {
-      //      low = min(low, bbi->lower());
-      //      upp = max(upp, bbi->upper());
-      //    }
+    // Loop over all levels finer than this one.
     
-      ivect low = bbi->lower();
-      ivect upp = bbi->upper();
-      
-      // low and upp now define the starting bbox.
+    leave_singlemap_mode(const_cast<cGH *> (cctkGH));
+    leave_level_mode(const_cast<cGH *> (cctkGH));
+    for (CCTK_INT rl = called_on_rl; rl < finest_current_rl; ++rl) {
+      enter_level_mode(const_cast<cGH *> (cctkGH), rl);
+      enter_singlemap_mode(const_cast<cGH *> (cctkGH), called_on_map);
 
-      ibbox bb(low, upp, bbs.at(0).stride());
-    
       if (verbose) {
         ostringstream buf;
-        buf << "Found the local size of the box: " << endl << bb;
+        buf << "Entering level " << rl << " map " << called_on_map;
         CCTK_INFO(buf.str().c_str());
       }
       
-      vector<CCTK_INT> mask(prod(bb.shape()/bb.stride()), 0);
+      // So the full algorithm should look something like:
       
-      if (veryverbose) {
-        ostringstream buf;
-        buf << "Allocated mask size: " << bb.shape()/bb.stride() 
-            << " (points: " << prod(bb.shape()/bb.stride()) << ")";
-        CCTK_INFO(buf.str().c_str());
-      }
+      // Find how big the first bounding box should be on this level
+      //   Do this by finding min lower and max upper bounds of all bboxes
+      // Allocate box
+      // Fill errors from local arrays
+      // If grandchildren exist use their bboxes (expanded) to add to errors
+      // Reduce errors (MPI sum)
+      // Set errors to 0/1
+      // Define vectors of bboxes final (empty) and todo (contains bbox)
+      // Define vector of masks (contains error mask)
+      // Loop over all entries in todo:
+      //   Setup appropriate 1d array memory
+      //   Call fortran routine
+      //   If return is:
+      //     zero: add bbox to final
+      //     one:  add new bbox to todo and assoc mask to masklist
+      //     two:  add both new bboxs to todo and assoc masks to masklist
+      //   
       
-
-      // Setup the mask.
+      vector<ibbox> bbs = local_bbsss.at(mglevel).at(reflevel);
       
-      const ibbox& baseext = 
-        vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior;
-      ivect imin = (bb.lower() - baseext.lower())/bb.stride(), 
-        imax = (bb.upper() - baseext.lower())/bb.stride();
-
-      BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) 
-      {
-        const CCTK_REAL *error_var_ptr = 
-          static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, error_var));
-        const CCTK_REAL *x_var_ptr = 
-          static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::x"));
-        const CCTK_REAL *y_var_ptr = 
-          static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::y"));
-        const CCTK_REAL *z_var_ptr = 
-          static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::z"));
-     
-        assert(all(imin >= 0));
-        assert(all(imax >= 0));
-        // FIXME: Why should the following assert be true?
-        //      assert(all(imax < ivect::ref(cctkGH->cctk_lsh)));
-        assert(all(imin <= imax));
+      stack<ibbox> final;
       
-        for (CCTK_INT k = 0; k < cctkGH->cctk_lsh[2]; ++k) {
-          for (CCTK_INT j = 0; j < cctkGH->cctk_lsh[1]; ++j) {
-            for (CCTK_INT i = 0; i < cctkGH->cctk_lsh[0]; ++i) {
-              CCTK_INT index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-              CCTK_REAL local_error = abs(error_var_ptr[index]);
-              if (local_error > max_error) {
-                CCTK_INT ii = i + cctkGH->cctk_lbnd[0] - imin[0];
-                CCTK_INT jj = j + cctkGH->cctk_lbnd[1] - imin[1];
-                CCTK_INT kk = k + cctkGH->cctk_lbnd[2] - imin[2];
-                // Check that this point actually intersects with 
-                // this box (if this component was actually a
-                // different grid on the same processor, it need not)
-                if ( (ii >= 0) and (jj >= 0) and (kk >= 0) and 
-                     (ii <= imax[0] - imin[0]) and
-                     (jj <= imax[1] - imin[1]) and
-                     (kk <= imax[2] - imin[2]) )
-                {
-                  assert (ii >= 0);
-                  assert (jj >= 0);
-                  assert (kk >= 0);
-                  assert (ii <= imax[0] - imin[0]);
-                  assert (jj <= imax[1] - imin[1]);
-                  assert (kk <= imax[2] - imin[2]);
-                  CCTK_INT mindex = ii + 
-                    (imax[0] - imin[0] + 1)*(jj + (imax[1] - imin[1] + 1) * kk);
-                  mask[mindex] = 1;
-                  if (veryverbose) {
-                    CCTK_VInfo(CCTK_THORNSTRING, "In error at point"
-                               "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]",
-                               x_var_ptr[index],
-                               y_var_ptr[index],
-                               z_var_ptr[index],
-                               ii, jj, kk, i,j,k);
-                  }
-                }
-              }
-            }
-          }
-        }
-      } END_LOCAL_COMPONENT_LOOP;
+      vector<vector<ibbox> > bbss = local_bbsss.at(0);
       
-      // FIXME    
-      // If grandchildren exist use their bboxes (expanded) to add to errors
-      // Reduce errors (MPI sum)
-      // Set errors to 0/1
-
-      // Pad the errors: stage 1 - buffer points marked as 2.
+      bool did_regrid = false;
+            
+      rvect physical_min, physical_max;
+      rvect interior_min, interior_max;
+      rvect exterior_min, exterior_max;
+      rvect base_spacing;
+      int ierr = GetDomainSpecification
+        (dim, &physical_min[0], &physical_max[0],
+         &interior_min[0], &interior_max[0],
+         &exterior_min[0], &exterior_max[0], &base_spacing[0]);
+      assert (!ierr);
+      
+      for ( vector<ibbox>::const_iterator bbi = bbs.begin(); 
+            bbi != bbs.end();
+            ++bbi) 
+      {
+        
+        ivect low = bbi->lower();
+        ivect upp = bbi->upper();
+        
+        // low and upp now define the starting bbox.
+        
+        ibbox bb(low, upp, bbs.at(0).stride());
+        
+        if (verbose) {
+          ostringstream buf;
+          buf << "Found the local size of the box: " << endl << bb;
+          CCTK_INFO(buf.str().c_str());
+        }
+        
+        vector<CCTK_INT> mask(prod(bb.shape()/bb.stride()), 0);
+        
+        if (veryverbose) {
+          ostringstream buf;
+          buf << "Allocated mask size: " << bb.shape()/bb.stride() 
+              << " (points: " << prod(bb.shape()/bb.stride()) << ")";
+          CCTK_INFO(buf.str().c_str());
+        }
+        
+        // Setup the mask.
+        
+        const ibbox& baseext = 
+          vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior;
+        ivect imin = (bb.lower() - baseext.lower())/bb.stride(), 
+          imax = (bb.upper() - baseext.lower())/bb.stride();
+        
+        BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) 
+        {
+          const CCTK_REAL *error_var_ptr = 
+            static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                          0, error_var));
+          const CCTK_REAL *x_var_ptr = 
+            static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                          0, "Grid::x"));
+          const CCTK_REAL *y_var_ptr = 
+            static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                          0, "Grid::y"));
+          const CCTK_REAL *z_var_ptr = 
+            static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                          0, "Grid::z"));
+          
 
-      for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
-        for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) {
-          for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) {
-            CCTK_INT index =  i + 
-              (imax[0] - imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k);
-            if (mask[index] == 1) {
-              for (CCTK_INT kk = max(k - pad, 0); 
-                   kk < min(k + pad + 1, imax[2] - imin[2] + 1);
-                   ++kk)
-              {
-                for (CCTK_INT jj = max(j - pad, 0); 
-                     jj < min(j + pad + 1, imax[1] - imin[1] + 1);
-                     ++jj)
-                {
-                  for (CCTK_INT ii = max(i - pad, 0); 
-                       ii < min(i + pad + 1, imax[0] - imin[0] + 1);
-                       ++ii)
+          // These can actually be negative if the parent shrinks.
+          // Of course, the final grid should still be properly
+          // nested...
+          //          assert(all(imin >= 0));
+          //          assert(all(imax >= 0));
+          // FIXME: Why should the following assert be true?
+          //      assert(all(imax < ivect::ref(cctkGH->cctk_lsh)));
+          assert(all(imin <= imax));
+          
+          for (CCTK_INT k = 0; k < cctkGH->cctk_lsh[2]; ++k) {
+            for (CCTK_INT j = 0; j < cctkGH->cctk_lsh[1]; ++j) {
+              for (CCTK_INT i = 0; i < cctkGH->cctk_lsh[0]; ++i) {
+                CCTK_INT index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL local_error = abs(error_var_ptr[index]);
+                if (local_error > max_error) {
+                  CCTK_INT ii = i + cctkGH->cctk_lbnd[0] - imin[0];
+                  CCTK_INT jj = j + cctkGH->cctk_lbnd[1] - imin[1];
+                  CCTK_INT kk = k + cctkGH->cctk_lbnd[2] - imin[2];
+                  // Check that this point actually intersects with 
+                  // this box (if this component was actually a
+                  // different grid on the same processor, it need not)
+                  if ( (ii >= 0) and (jj >= 0) and (kk >= 0) and 
+                       (ii <= imax[0] - imin[0]) and
+                       (jj <= imax[1] - imin[1]) and
+                       (kk <= imax[2] - imin[2]) )
                   {
+                    assert (ii >= 0);
+                    assert (jj >= 0);
+                    assert (kk >= 0);
+                    assert (ii <= imax[0] - imin[0]);
+                    assert (jj <= imax[1] - imin[1]);
+                    assert (kk <= imax[2] - imin[2]);
                     CCTK_INT mindex = ii + 
                       (imax[0] - imin[0] + 1)*
                       (jj + (imax[1] - imin[1] + 1) * kk);
-                    if (!mask[mindex]) mask[mindex] = 2;
+                    mask[mindex] = 1;
+                    if (veryverbose) {
+                      CCTK_VInfo(CCTK_THORNSTRING, "In error at point"
+                                 "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]",
+                                 x_var_ptr[index],
+                                 y_var_ptr[index],
+                                 z_var_ptr[index],
+                                 ii, jj, kk, i,j,k);
+                    }
                   }
                 }
               }
             }
           }
-        }
-      }
-      // stage 2: all buffer points marked truly in error.
-      // Also mark if there are any errors.
-      bool should_regrid = false;
-      for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
-        for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) {
-          for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) {
-            CCTK_INT index =  i + 
-              (imax[0]-imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k);
-            if (mask[index] > 1) mask[index] = 1;
-            if ((veryverbose) and (mask[index])) {
-              CCTK_VInfo(CCTK_THORNSTRING, "Mask set at point"
-                         "\n[%d,%d,%d]",
-                         i,j,k);
-            } 
-            should_regrid |= (mask[index]);
-            did_regrid |= should_regrid;
+        } END_LOCAL_COMPONENT_LOOP;
+        
+        // Instead check the error on child level, if exists
+        // This should fix the "orphaned grandchild" problem
+        
+        if (bbss.size() > reflevel+1) {
+                
+          CCTK_INT currentml = mglevel;
+          CCTK_INT currentrl = reflevel;
+          CCTK_INT currentmap = carpetGH.map;            
+          
+          leave_singlemap_mode(const_cast<cGH *> (cctkGH));
+          leave_level_mode(const_cast<cGH *> (cctkGH));
+      
+          enter_level_mode(const_cast<cGH *> (cctkGH), currentrl + 1);
+          enter_singlemap_mode(const_cast<cGH *> (cctkGH), currentmap);
+                          
+          if (verbose) {
+            ostringstream buf;
+            buf << "Checking for errors on child level " 
+                << reflevel << " map " << currentmap;
+            CCTK_INFO(buf.str().c_str());
           }
-        }
-      }    
       
-      // Define vectors of bboxes final (empty) and todo (contains bbox)
+          BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+            const CCTK_REAL *error_var_ptr = 
+              static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                            0, error_var));
+            const CCTK_REAL *x_var_ptr = 
+              static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                            0, "Grid::x"));
+            const CCTK_REAL *y_var_ptr = 
+              static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                            0, "Grid::y"));
+            const CCTK_REAL *z_var_ptr = 
+              static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 
+                                                            0, "Grid::z"));
+            
+            //            assert(all(imin >= 0));
+            //            assert(all(imax >= 0));
+            // FIXME: Why should the following assert be true?
+            //      assert(all(imax < ivect::ref(cctkGH->cctk_lsh)));
+            assert(all(imin <= imax));
+            
+            for (CCTK_INT k = 0; k < cctkGH->cctk_lsh[2]; ++k) {
+              for (CCTK_INT j = 0; j < cctkGH->cctk_lsh[1]; ++j) {
+                for (CCTK_INT i = 0; i < cctkGH->cctk_lsh[0]; ++i) {
+                  CCTK_INT index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                  CCTK_REAL local_error = abs(error_var_ptr[index]);
+                  if (local_error > max_error) {
+                    CCTK_INT ii = i + cctkGH->cctk_lbnd[0] - imin[0];
+                    CCTK_INT jj = j + cctkGH->cctk_lbnd[1] - imin[1];
+                    CCTK_INT kk = k + cctkGH->cctk_lbnd[2] - imin[2];
+                    // Check that this point actually intersects with 
+                    // this box (if this component was actually a
+                    // different grid on the same processor, it need not)
+                    if ( (ii >= 0) and (jj >= 0) and (kk >= 0) and 
+                         (ii <= imax[0] - imin[0]) and
+                         (jj <= imax[1] - imin[1]) and
+                         (kk <= imax[2] - imin[2]) )
+                    {
+                      assert (ii >= 0);
+                      assert (jj >= 0);
+                      assert (kk >= 0);
+                      assert (ii <= imax[0] - imin[0]);
+                      assert (jj <= imax[1] - imin[1]);
+                      assert (kk <= imax[2] - imin[2]);
+                      CCTK_INT mindex = ii + 
+                        (imax[0] - imin[0] + 1)*
+                        (jj + (imax[1] - imin[1] + 1) * kk);
+                      mask[mindex] = 1;
+                      if (veryverbose) {
+                        CCTK_VInfo(CCTK_THORNSTRING, "In error at point"
+                                   "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]",
+                                   x_var_ptr[index],
+                                   y_var_ptr[index],
+                                   z_var_ptr[index],
+                                   ii, jj, kk, i,j,k);
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          } END_LOCAL_COMPONENT_LOOP;
+        
+          leave_singlemap_mode(const_cast<cGH *> (cctkGH));
+          leave_level_mode(const_cast<cGH *> (cctkGH));
+          
+          enter_level_mode(const_cast<cGH *> (cctkGH), currentrl);
+          enter_singlemap_mode(const_cast<cGH *> (cctkGH), currentmap);
 
-      if (should_regrid) {
-    
-        stack<ibbox> todo;
+        }
+        
+        // FIXME
         
-        todo.push(bb);
+        // Reduce errors (MPI sum)
+        // Set errors to 0/1
         
-        // Define vector of masks (contains error mask)
+        // Pad the errors: stage 1 - buffer points marked as 2.
         
-        stack<vector<CCTK_INT> > masklist;
+        for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
+          for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) {
+            for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) {
+              CCTK_INT index =  i + 
+                (imax[0] - imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k);
+              if (mask[index] == 1) {
+                for (CCTK_INT kk = max(k - pad, 0); 
+                     kk < min(k + pad + 1, imax[2] - imin[2] + 1);
+                       ++kk)
+                {
+                  for (CCTK_INT jj = max(j - pad, 0); 
+                       jj < min(j + pad + 1, imax[1] - imin[1] + 1);
+                       ++jj)
+                  {
+                    for (CCTK_INT ii = max(i - pad, 0); 
+                         ii < min(i + pad + 1, imax[0] - imin[0] + 1);
+                         ++ii)
+                    {
+                      CCTK_INT mindex = ii + 
+                        (imax[0] - imin[0] + 1)*
+                        (jj + (imax[1] - imin[1] + 1) * kk);
+                      if (!mask[mindex]) mask[mindex] = 2;
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+        // stage 2: all buffer points marked truly in error.
+        // Also mark if there are any errors.
+        bool should_regrid = false;
+        for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
+          for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) {
+            for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) {
+              CCTK_INT index =  i + 
+                (imax[0]-imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k);
+              if (mask[index] > 1) mask[index] = 1;
+              if ((veryverbose) and (mask[index])) {
+                CCTK_VInfo(CCTK_THORNSTRING, "Mask set at point"
+                           "\n[%d,%d,%d]",
+                           i,j,k);
+              } 
+              should_regrid |= (mask[index]);
+              did_regrid |= should_regrid;
+            }
+          }
+        }    
         
-        masklist.push(mask);
+        if (verbose) {
+          ostringstream buf;
+          buf << "Finished looking for errors on level " 
+              << reflevel << endl << "should_regrid " << should_regrid
+              << " did_regrid " << did_regrid;
+          CCTK_INFO(buf.str().c_str());
+        }
         
-        // Loop over all entries in todo:
-        //   Setup appropriate 1d array memory
-        //   Call fortran routine
-        //   If return is:
-        //     zero: add bbox to final
-        //     one:  add new bbox to todo and assoc mask to masklist
-        //     two:  add both new bboxs to todo and assoc masks to masklist
+        // Define vectors of bboxes final (empty) and todo (contains bbox)
         
-        while (!todo.empty())
-        {
-          
-          ibbox bb = todo.top(); todo.pop();
-          vector<CCTK_INT> mask = masklist.top(); masklist.pop();
-          
-          CCTK_INT nx = bb.shape()[0]/bb.stride()[0];
-          CCTK_INT ny = bb.shape()[1]/bb.stride()[1];
-          CCTK_INT nz = bb.shape()[2]/bb.stride()[2];
-          
-          if (verbose) {
-            ostringstream buf;
-            buf << "todo loop. Box: " << endl << bb;
-            CCTK_INFO(buf.str().c_str());
-          }
+        if (should_regrid) {
           
-          vector<CCTK_INT> sum_x(nx, 0);
-          vector<CCTK_INT> sig_x(nx, 0);
-          vector<CCTK_INT> sum_y(ny, 0);
-          vector<CCTK_INT> sig_y(ny, 0);
-          vector<CCTK_INT> sum_z(nz, 0);
-          vector<CCTK_INT> sig_z(nz, 0);
+          stack<ibbox> todo;
           
-          CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3];
+          todo.push(bb);
           
-          for (CCTK_INT d = 0; d < 3; ++d) {
-            fbbox[0][d] = bb.lower()[d];
-            fbbox[1][d] = bb.upper()[d];
-            fbbox[2][d] = bb.stride()[d];
-          }
+          // Define vector of masks (contains error mask)
           
-          CCTK_INT didit;
+          stack<vector<CCTK_INT> > masklist;
           
-          CCTK_FNAME(check_box)(nx, ny, nz,
-                                &mask.front(),
-                                &sum_x.front(), &sum_y.front(), &sum_z.front(),
-                                &sig_x.front(), &sig_y.front(), &sig_z.front(),
-                                fbbox,
-                                fbbox1, fbbox2,
-                                min_width, min_fraction,
-                                didit);
+          masklist.push(mask);
           
-          if (didit == 0) {
+          // Loop over all entries in todo:
+          //   Setup appropriate 1d array memory
+          //   Call fortran routine
+          //   If return is:
+          //     zero: add bbox to final
+          //     one:  add new bbox to todo and assoc mask to masklist
+          //     two:  add both new bboxs to todo and assoc masks to masklist
           
-            final.push(bb);          
-
+          while (!todo.empty())
+          {
+            
+            ibbox bb = todo.top(); todo.pop();
+            vector<CCTK_INT> mask = masklist.top(); masklist.pop();
+            
+            CCTK_INT nx = bb.shape()[0]/bb.stride()[0];
+            CCTK_INT ny = bb.shape()[1]/bb.stride()[1];
+            CCTK_INT nz = bb.shape()[2]/bb.stride()[2];
+            
             if (verbose) {
               ostringstream buf;
-              buf << "todo loop. Box pushed to final: " 
-                  << endl << bb;
+              buf << "todo loop. Box: " << endl << bb;
               CCTK_INFO(buf.str().c_str());
             }
-          }
-          else if (didit == 1) {
             
-            ibbox newbbox1(ivect::ref(&fbbox1[0][0]),
-                           ivect::ref(&fbbox1[1][0]),
-                           ivect::ref(&fbbox1[2][0]));
-            todo.push(newbbox1);
+            vector<CCTK_INT> sum_x(nx, 0);
+            vector<CCTK_INT> sig_x(nx, 0);
+            vector<CCTK_INT> sum_y(ny, 0);
+            vector<CCTK_INT> sig_y(ny, 0);
+            vector<CCTK_INT> sum_z(nz, 0);
+            vector<CCTK_INT> sig_z(nz, 0);
             
-            CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0];
-            CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1];
-            CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2];
+            CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3];
             
-            vector<CCTK_INT>  
-              newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0);
-            
-            CCTK_FNAME(copy_mask)(nx, ny, nz,
-                                  &mask.front(), fbbox,
-                                  dnx, dny, dnz,
-                                  &newmask1.front(), fbbox1);
-            masklist.push(newmask1);
+            for (CCTK_INT d = 0; d < 3; ++d) {
+              fbbox[0][d] = bb.lower()[d];
+              fbbox[1][d] = bb.upper()[d];
+              fbbox[2][d] = bb.stride()[d];
+            }
+              
+            CCTK_INT didit;
+              
+            CCTK_FNAME(check_box)(nx, ny, nz,
+                                  &mask.front(),
+                                  &sum_x.front(), &sum_y.front(),
+                                  &sum_z.front(),
+                                  &sig_x.front(), &sig_y.front(),
+                                  &sig_z.front(),
+                                  fbbox,
+                                  fbbox1, fbbox2,
+                                  min_width, min_fraction,
+                                  didit);
             
-            if (verbose) {
-              ostringstream buf;
-              buf << "todo loop. New (single) box created: " 
-                  << endl << newbbox1;
-              CCTK_INFO(buf.str().c_str());
+            if (didit == 0) {
+              
+              final.push(bb);          
+              
+              if (verbose) {
+                ostringstream buf;
+                buf << "todo loop. Box pushed to final: " 
+                    << endl << bb;
+                CCTK_INFO(buf.str().c_str());
+              }
+            }
+            else if (didit == 1) {
+              
+              ibbox newbbox1(ivect::ref(&fbbox1[0][0]),
+                             ivect::ref(&fbbox1[1][0]),
+                             ivect::ref(&fbbox1[2][0]));
+              todo.push(newbbox1);
+              
+              CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0];
+              CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1];
+              CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2];
+              
+              vector<CCTK_INT>  
+                newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0);
+              
+              CCTK_FNAME(copy_mask)(nx, ny, nz,
+                                    &mask.front(), fbbox,
+                                    dnx, dny, dnz,
+                                    &newmask1.front(), fbbox1);
+              masklist.push(newmask1);
+              
+              if (verbose) {
+                ostringstream buf;
+                buf << "todo loop. New (single) box created: " 
+                    << endl << newbbox1;
+                CCTK_INFO(buf.str().c_str());
+              }
+            }
+            else if (didit == 2) {
+              
+              ibbox newbbox1(ivect::ref(&fbbox1[0][0]),
+                             ivect::ref(&fbbox1[1][0]),
+                             ivect::ref(&fbbox1[2][0]));
+              todo.push(newbbox1);
+              ibbox newbbox2(ivect::ref(&fbbox2[0][0]),
+                             ivect::ref(&fbbox2[1][0]),
+                             ivect::ref(&fbbox2[2][0]));
+              todo.push(newbbox2);
+              
+              CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0];
+              CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1];
+              CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2];
+              
+              vector<CCTK_INT>  
+                newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0);
+              
+              CCTK_FNAME(copy_mask)(nx, ny, nz,
+                                    &mask.front(), fbbox,
+                                    dnx, dny, dnz,
+                                    &newmask1.front(), fbbox1);
+              masklist.push(newmask1);
+              
+              dnx = newbbox2.shape()[0]/newbbox2.stride()[0];
+              dny = newbbox2.shape()[1]/newbbox2.stride()[1];
+              dnz = newbbox2.shape()[2]/newbbox2.stride()[2];
+              
+              vector<CCTK_INT>  
+                newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0);
+              
+              CCTK_FNAME(copy_mask)(nx, ny, nz,
+                                    &mask.front(), fbbox,
+                                    dnx, dny, dnz,
+                                    &newmask2.front(), fbbox2);
+              masklist.push(newmask2);
+              
+              if (verbose) {
+                ostringstream buf;
+                buf << "todo loop. New (double) box created. Box 1: " 
+                    << endl << newbbox1
+                    << "                                     Box 2: "
+                    << endl << newbbox2;
+                CCTK_INFO(buf.str().c_str());
+              }
+            }
+            else {
+              CCTK_WARN(0, "The fortran routine must be confused.");
             }
-          }
-          else if (didit == 2) {
             
-            ibbox newbbox1(ivect::ref(&fbbox1[0][0]),
-                           ivect::ref(&fbbox1[1][0]),
-                           ivect::ref(&fbbox1[2][0]));
-            todo.push(newbbox1);
-            ibbox newbbox2(ivect::ref(&fbbox2[0][0]),
-                           ivect::ref(&fbbox2[1][0]),
-                           ivect::ref(&fbbox2[2][0]));
-            todo.push(newbbox2);
+          } // loop over todo vector (boxes needing to be done).
+        } // should regrid.
+      } // Loop over boxes on the parent grid.
+      
+      if (did_regrid) {
+          
+        // Fixup the stride
+        vector<ibbox> newbbs;
+        vector<bbvect> obs;
+        while (! final.empty()) {
+          ibbox bb = final.top(); final.pop();
+          
+          if (veryverbose) {
+            ostringstream buf;
+            buf << "Looping over the final list. Box is:"
+                << endl << bb;
+            CCTK_INFO(buf.str().c_str());
+          }
+          
+          ivect ilo = bb.lower();
+          ivect ihi = bb.upper();
+          rvect lo = int2pos(cctkGH, hh, ilo, reflevel);
+          rvect hi = int2pos(cctkGH, hh, ihi, reflevel);
+          rvect str = base_spacing * 
+            ipow((CCTK_REAL)mgfact, basemglevel) /
+            ipow(reffact, reflevel);
+          rbbox newbbcoord(lo, hi, str);
+          
+          if (veryverbose) {
+            ostringstream buf;
+            buf << "Dealing with boundaries. Coord box is:"
+                << endl << newbbcoord;
+            CCTK_INFO(buf.str().c_str());
+          }
+          
+          // Set the correct ob here.
+          
+          bbvect ob(false);
+          for (int d=0; d<dim; ++d) {
+            assert (mglevel==0);
             
-            CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0];
-            CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1];
-            CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2];
+            // Find the size of the physical domain
             
-            vector<CCTK_INT>  
-              newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0);
+            rvect const spacing = base_spacing * 
+              ipow((CCTK_REAL)mgfact, basemglevel) /
+              ipow(reffact, reflevel+1);
+            ierr = ConvertFromPhysicalBoundary
+              (dim, &physical_min[0], &physical_max[0],
+               &interior_min[0], &interior_max[0],
+               &exterior_min[0], &exterior_max[0], &spacing[0]);
+            assert (!ierr);
             
-            CCTK_FNAME(copy_mask)(nx, ny, nz,
-                                  &mask.front(), fbbox,
-                                  dnx, dny, dnz,
-                                  &newmask1.front(), fbbox1);
-            masklist.push(newmask1);
+            // If need be clip the domain
             
-            dnx = newbbox2.shape()[0]/newbbox2.stride()[0];
-            dny = newbbox2.shape()[1]/newbbox2.stride()[1];
-            dnz = newbbox2.shape()[2]/newbbox2.stride()[2];
+            rvect lo = newbbcoord.lower();
+            if (newbbcoord.lower()[d] < physical_min[d]) {
+              lo[d] = exterior_min[d];
+            }
+            rvect up = newbbcoord.upper();
+            if (newbbcoord.upper()[d] > physical_max[d]) {
+              up[d] = exterior_max[d];
+            }
+            rvect str = newbbcoord.stride();
             
-            vector<CCTK_INT>  
-              newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0);
+            // Set the ob if outside the physical domain
             
-            CCTK_FNAME(copy_mask)(nx, ny, nz,
-                                  &mask.front(), fbbox,
-                                  dnx, dny, dnz,
-                                  &newmask2.front(), fbbox2);
-            masklist.push(newmask2);
+            ob[d][0] = 
+              abs(lo[d] - exterior_min[d]) < 1.0e-6 * spacing[d];
+            ob[d][1] = 
+              abs(up[d] - exterior_max[d]) < 1.0e-6 * spacing[d];
             
-            if (verbose) {
+            if (veryverbose) {
               ostringstream buf;
-              buf << "todo loop. New (double) box created. Box 1: " 
-                  << endl << newbbox1
-                  << "                                     Box 2: "
-                  << endl << newbbox2;
+              buf << "Done clipping domain:"
+                  << endl << lo << endl << up << endl << str;
               CCTK_INFO(buf.str().c_str());
+            } 
+            
+            // Check that the striding is correct.
+            
+            rvect remainder = (up - lo) / str - floor( (up - lo) / str );
+            
+            for (int d=0; d < dim; ++d) {
+              if ( abs(remainder[d]) > 1.e-6 ) {
+                if (ob[d][0]) {
+                  up[d] += str[d] * (1 - remainder[d]);
+                }
+                else if (ob[d][1]) {
+                  lo[d] -= str[d] * remainder[d];
+                }
+              }
             }
+            
+            if (veryverbose) {
+              ostringstream buf;
+              buf << "Corrected coords for striding:"
+                  << endl << lo << endl << up << endl << str;
+              CCTK_INFO(buf.str().c_str());
+            } 
+            
+            newbbcoord = rbbox(lo, up, str);
           }
-          else {
-            CCTK_WARN(0, "The fortran routine must be confused.");
-          }
-          
-        } // loop over todo vector (boxes needing to be done).
-      } // should regrid.
-    } // Loop over boxes on the parent grid.
-
-    if (did_regrid) {
-
-      // Fixup the stride
-      vector<ibbox> newbbs;
-      vector<bbvect> obs;
-      while (! final.empty()) {
-        ibbox bb = final.top(); final.pop();
-        
-        if (veryverbose) {
-          ostringstream buf;
-          buf << "Looping over the final list. Box is:"
-              << endl << bb;
-          CCTK_INFO(buf.str().c_str());
-        }
-
-        ivect ilo = bb.lower();
-        ivect ihi = bb.upper();
-        rvect lo = int2pos(cctkGH, hh, ilo, reflevel);
-        rvect hi = int2pos(cctkGH, hh, ihi, reflevel);
-        rvect str = base_spacing * 
-          ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel);
-        rbbox newbbcoord(lo, hi, str);
-        
-        if (veryverbose) {
-          ostringstream buf;
-          buf << "Dealing with boundaries. Coord box is:"
-              << endl << newbbcoord;
-          CCTK_INFO(buf.str().c_str());
-        }
-        
-        // FIXME: Set the correct ob here.
-        
-        bbvect ob(false);
-        for (int d=0; d<dim; ++d) {
-          assert (mglevel==0);
-          
-          // Find the size of the physical domain
-          
-          rvect const spacing = base_spacing * 
-            ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel+1);
-          ierr = ConvertFromPhysicalBoundary
-            (dim, &physical_min[0], &physical_max[0],
-             &interior_min[0], &interior_max[0],
-             &exterior_min[0], &exterior_max[0], &spacing[0]);
-          assert (!ierr);
-          
-          // If need be clip the domain
-          
-          rvect lo = newbbcoord.lower();
-          if (newbbcoord.lower()[d] < physical_min[d]) {
-            lo[d] = exterior_min[d];
-          }
-          rvect up = newbbcoord.upper();
-          if (newbbcoord.upper()[d] > physical_max[d]) {
-            up[d] = exterior_max[d];
-          }
-          rvect str = newbbcoord.stride();
-          
-          // Set the ob if outside the physical domain
-          
-          ob[d][0] = 
-            abs(lo[d] - exterior_min[d]) < 1.0e-6 * spacing[d];
-          ob[d][1] = 
-            abs(up[d] - exterior_max[d]) < 1.0e-6 * spacing[d];
-
-          if (veryverbose) {
+          if (verbose) {
             ostringstream buf;
-            buf << "Done clipping domain:"
-                << endl << lo << endl << up << endl << str;
+            buf << "Done dealing with boundaries. Coord box is:"
+                << endl << newbbcoord << endl
+                << "obox is:" << endl << ob;
             CCTK_INFO(buf.str().c_str());
-          } 
+          }
           
-          // Check that the striding is correct.
+          // Convert back to integer coordinates
+          // We have to do this on the fine grid to ensure that
+          // it is correct for an outer boundary with odd numbers
+          // of ghost zones where the bbox does not align with the parent.
           
-          rvect remainder = (up - lo) / str - floor( (up - lo) / str );
-
+          ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel+1);
+          ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel+1);
+          ivect istr = bb.stride() / reffact;
+          
+          // Check that the width is sufficient
+          // This can only be too small if the domain was clipped
           for (int d=0; d < dim; ++d) {
-            if ( abs(remainder[d]) > 1.e-6 ) {
+            if (ihi[d] - ilo[d] < min_width * istr[d]) {
               if (ob[d][0]) {
-                up[d] += str[d] * (1 - remainder[d]);
+                if (ob[d][1]) {
+                  CCTK_WARN(0, "The domain is too small?!");
+                }
+                ihi[d] = ilo[d] + min_width * istr[d];
               }
               else if (ob[d][1]) {
-                lo[d] -= str[d] * remainder[d];
+                if (ob[d][0]) {
+                  CCTK_WARN(0, "The domain is too small?!");
+                }
+                ilo[d] = ihi[d] - min_width * istr[d];
+              }
+              else {
+                ostringstream buf;
+                buf << "The grid is unclipped and too small?" << endl 
+                    << ilo << endl << ihi << endl << istr << endl << d;
+                CCTK_WARN(0, buf.str().c_str());
               }
             }
           }
           
           if (veryverbose) {
             ostringstream buf;
-            buf << "Corrected coords for striding:"
-                << endl << lo << endl << up << endl << str;
+            buf << "Corrected integer coords for min_width:"
+                << endl << ilo << endl << ihi << endl << istr;
             CCTK_INFO(buf.str().c_str());
           } 
-
-          newbbcoord = rbbox(lo, up, str);
-        }
-        if (verbose) {
-          ostringstream buf;
-          buf << "Done dealing with boundaries. Coord box is:"
-              << endl << newbbcoord << endl
-              << "obox is:" << endl << ob;
-          CCTK_INFO(buf.str().c_str());
+          
+          ibbox newbb(ilo, ihi, istr);          
+          
+          if (verbose) {
+            ostringstream buf;
+            buf << "After dealing with boundaries. Final box is:"
+                << endl << newbb;
+            CCTK_INFO(buf.str().c_str());
+          }
+          
+          newbbs.push_back (newbb);
+          obs.push_back(ob);
         }
         
-        // Convert back to integer coordinates
-        // We have to do this on the fine grid to ensure that
-        // it is correct for an outer boundary with odd numbers
-        // of ghost zones where the bbox does not align with the parent.
         
-        ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel+1);
-        ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel+1);
-        ivect istr = bb.stride() / reffact;
+        // FIXME: check if the newbbs is really different
+        // from the current bbs
+        //        if not, set do_recompose = 0
+        bbs = newbbs;
         
-        // Check that the width is sufficient
-        // This can only be too small if the domain was clipped
-        for (int d=0; d < dim; ++d) {
-          if (ihi[d] - ilo[d] < min_width * istr[d]) {
-            if (ob[d][0]) {
-              if (ob[d][1]) {
-                CCTK_WARN(0, "The domain is too small?!");
-              }
-              ihi[d] = ilo[d] + min_width * istr[d];
-            }
-            else if (ob[d][1]) {
-              if (ob[d][0]) {
-                CCTK_WARN(0, "The domain is too small?!");
-              }
-              ilo[d] = ihi[d] - min_width * istr[d];
-            }
-            else {
-              CCTK_WARN(0, "The grid is unclipped and too small?!");
-            }
+        // Set local bbss
+        
+        if (bbss.size() < reflevel+2) {
+          if (verbose) {
+            CCTK_INFO("Adding new refinement level");
           }
+          bbss.resize(reflevel+2);
+          obss.resize(reflevel+2);
+          pss.resize(reflevel+2);
         }
-
-        if (veryverbose) {
-          ostringstream buf;
-          buf << "Corrected integer coords for min_width:"
-              << endl << ilo << endl << ihi << endl << istr;
-          CCTK_INFO(buf.str().c_str());
-        } 
+        bbss.at(reflevel+1) = bbs;
+        obss.at(reflevel+1) = obs;
+        MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
         
-        ibbox newbb(ilo, ihi, istr);          
+        // make multiprocessor aware
+        gh::cprocs ps;
+        SplitRegions (cctkGH, bbs, obs, ps);    
         
-        if (verbose) {
-          ostringstream buf;
-          buf << "After dealing with boundaries. Final box is:"
-              << endl << newbb;
-          CCTK_INFO(buf.str().c_str());
-        }
+        bbss.at(reflevel+1) = bbs;
+        obss.at(reflevel+1) = obs;
+        pss.at(reflevel+1) = ps;
         
-        newbbs.push_back (newbb);
-        obs.push_back(ob);
-      }
-      
-      
-      // FIXME: check if the newbbs is really different from the current bbs
-      //        if not, set do_recompose = 0
-      bbs = newbbs;
-      
-      // Set local bbss
-      
-      if (bbss.size() < reflevel+2) {
-        if (verbose) {
-          CCTK_INFO("Adding new refinement level");
-        }
-        bbss.resize(reflevel+2);
-        obss.resize(reflevel+2);
-        pss.resize(reflevel+2);
-      }
-      bbss.at(reflevel+1) = bbs;
-      obss.at(reflevel+1) = obs;
-      MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
-      
-      // make multiprocessor aware
-      gh::cprocs ps;
-      SplitRegions (cctkGH, bbs, obs, ps);    
-      
-      bbss.at(reflevel+1) = bbs;
-      obss.at(reflevel+1) = obs;
-      pss.at(reflevel+1) = ps;
-    
-    } // did_regrid?
-    else
-    {
-      if (bbss.size() > reflevel+1) {
-        if (verbose) {
-          CCTK_INFO("Removing refinement level");
+      } // did_regrid?
+      else
+      {
+        if (bbss.size() > reflevel+1) {
+          if (verbose) {
+            CCTK_INFO("Removing refinement level");
+          }
         }
+        bbss.resize(reflevel+1);
+        obss.resize(reflevel+1);
+        // Set local bbsss
+        MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+        
+        pss.resize(reflevel+1);
+        
+        do_recompose = 1;
       }
-      bbss.resize(reflevel+1);
-      obss.resize(reflevel+1);
-      // Set local bbsss
-      MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+        
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
       
-      pss.resize(reflevel+1);
+      leave_singlemap_mode(const_cast<cGH *> (cctkGH));
+      leave_level_mode(const_cast<cGH *> (cctkGH));
       
-      do_recompose = 1;
-    }
-    
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+    } 
+
+    enter_level_mode(const_cast<cGH *> (cctkGH), called_on_rl);
+    enter_singlemap_mode(const_cast<cGH *> (cctkGH), called_on_map);
 
-//     cout << "bbsss done:" << endl << bbsss << endl;
-//     cout << "obss done:" << endl << obss << endl;
-//     cout << "pss done:" << endl << pss << endl;
+    if (verbose) {
+      ostringstream buf;
+      buf << "Done with it all. Total list is:"
+          << endl << local_bbsss;
+      CCTK_INFO(buf.str().c_str());
+    }      
     
     return do_recompose;
   }