Multilevel AMR

Basic support for more than 2 refinement levels with AMR. Basic tests
done with shock tubes seem to show that it works. Note that proper
nesting is not ensured as yet as we don't check for orphaned
grandchildren.

darcs-hash:20050218120448-58c7f-5d2e6605ac05f1363596c0e61f6026e8fc057ff8.gz

1 files changed, 325 insertions, 284 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index ce59cb533..eb18a5ca6 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -28,6 +28,8 @@ namespace CarpetAdaptiveRegrid {
   using namespace std;
   using namespace Carpet;
 
+  static gh::mexts local_bbsss;
+
   extern "C" {
     void CCTK_FCALL CCTK_FNAME(copy_mask)
       (const CCTK_INT& snx, const CCTK_INT& sny, const CCTK_INT& snz,
@@ -67,6 +69,15 @@ namespace CarpetAdaptiveRegrid {
     gh const & hh = *vhh.at(Carpet::map);
     
     assert (is_singlemap_mode());
+
+    if (local_bbsss.empty()) { // It's the first call
+      // Is this really the right thing to do on 
+      // multiprocessors?
+      local_bbsss = bbsss;
+    }
+
+    // FIXME: We should check that the local reflevel "agrees"
+    // with what is passed in.
     
     // In force mode (force == true) we do not check the
     // CarpetAdaptiveregrid parameters
@@ -129,44 +140,15 @@ namespace CarpetAdaptiveRegrid {
     //     two:  add both new bboxs to todo and assoc masks to masklist
     //   
 
-    vector<ibbox> bbs = bbsss.at(mglevel).at(reflevel);
+    vector<ibbox> bbs = local_bbsss.at(mglevel).at(reflevel);
     
-    ivect low = 100000, upp = -1;
-
-    for ( vector<ibbox>::const_iterator bbi = bbs.begin(); 
-          bbi != bbs.end();
-          ++bbi) 
-    {
-      low = min(low, bbi->lower());
-      upp = max(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());
-    }
-
+    stack<ibbox> final;
 
-    // Setup the mask.
+    vector<vector<ibbox> > bbss = local_bbsss.at(0);
       
-    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();
+    //    ivect low = 100000, upp = -1;
+
+    bool did_regrid = false;
 
     rvect physical_min, physical_max;
     rvect interior_min, interior_max;
@@ -178,268 +160,315 @@ namespace CarpetAdaptiveRegrid {
        &exterior_min[0], &exterior_max[0], &base_spacing[0]);
     assert (!ierr);
 
-    BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) 
+    for ( vector<ibbox>::const_iterator bbi = bbs.begin(); 
+          bbi != bbs.end();
+          ++bbi) 
     {
-      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"));
+      //      low = min(low, bbi->lower());
+      //      upp = max(upp, bbi->upper());
+      //    }
+    
+      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"));
      
-      assert(all(imin >= 0));
-      assert(all(imax >= 0));
-      assert(all(imin < ivect::ref(cctkGH->cctk_lsh)));
-      assert(all(imax < ivect::ref(cctkGH->cctk_lsh)));
-      assert(all(imin <= imax));
+        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);
-            if (abs(error_var_ptr[index]) > 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];
-              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);
+        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);
+              if (abs(error_var_ptr[index]) > 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;
-
-    // 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.
-
-    for (CCTK_INT k = 0; k < imax[2] - imin[2]; k++) {
-      for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) {
-        for (CCTK_INT i = 0; i < imax[0] - imin[0]; 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]);
-                 ++kk)
-            {
-              for (CCTK_INT jj = max(j-pad, 0); 
-                   jj < min(j+pad+1, imax[1] - imin[1]);
-                   ++jj)
+      } END_LOCAL_COMPONENT_LOOP;
+      
+      // 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.
+
+      for (CCTK_INT k = 0; k < imax[2] - imin[2]; k++) {
+        for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) {
+          for (CCTK_INT i = 0; i < imax[0] - imin[0]; 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]);
+                   ++kk)
               {
-                for (CCTK_INT ii = max(i-pad, 0); 
-                     ii < min(i+pad+1, imax[0] - imin[0]);
-                     ++ii)
+                for (CCTK_INT jj = max(j-pad, 0); 
+                     jj < min(j+pad+1, imax[1] - imin[1]);
+                     ++jj)
                 {
-                  CCTK_INT mindex = ii + 
-                    (imax[0] - imin[0] + 1)*
-                    (jj + (imax[1] - imin[1] + 1) * kk);
-                  if (!mask[mindex]) mask[mindex] = 2;
+                  for (CCTK_INT ii = max(i-pad, 0); 
+                       ii < min(i+pad+1, imax[0] - imin[0]);
+                       ++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]; k++) {
-      for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) {
-        for (CCTK_INT i = 0; i < imax[0] - imin[0]; 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]);
+      // 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]; k++) {
+        for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) {
+          for (CCTK_INT i = 0; i < imax[0] - imin[0]; 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;
+          }
         }
-      }
-    }    
-
-    // Define vectors of bboxes final (empty) and todo (contains bbox)
-      
-    stack<ibbox> final;
-
-    vector<vector<ibbox> > bbss = bbsss.at(0);
+      }    
       
+      // Define vectors of bboxes final (empty) and todo (contains bbox)
 
-    if (should_regrid) {
+      if (should_regrid) {
     
-      stack<ibbox> todo;
-  
-      todo.push(bb);
-      
-      // Define vector of masks (contains error mask)
-      
-      stack<vector<CCTK_INT> > masklist;
-      
-      masklist.push(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
-      
-      while (!todo.empty())
-      {
-      
-        ibbox bb = todo.top(); todo.pop();
-        vector<CCTK_INT> mask = masklist.top(); masklist.pop();
+        stack<ibbox> todo;
         
-        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];
+        todo.push(bb);
         
-        if (verbose) {
-          ostringstream buf;
-          buf << "todo loop. Box: " << endl << bb;
-          CCTK_INFO(buf.str().c_str());
-        }
-
-        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);
+        // Define vector of masks (contains error mask)
         
-        CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3];
+        stack<vector<CCTK_INT> > masklist;
         
-        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];
-        }
+        masklist.push(mask);
         
-        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);
+        // 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
         
-        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) {
+        while (!todo.empty())
+        {
           
-          ibbox newbbox1(ivect::ref(&fbbox1[0][0]),
-                         ivect::ref(&fbbox1[1][0]),
-                         ivect::ref(&fbbox1[2][0]));
-          todo.push(newbbox1);
+          ibbox bb = todo.top(); todo.pop();
+          vector<CCTK_INT> mask = masklist.top(); masklist.pop();
           
-          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 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];
           
-          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;
+            buf << "todo loop. Box: " << endl << bb;
             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);
+          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);
+          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_FNAME(copy_mask)(nx, ny, nz,
-                                &mask.front(), fbbox,
-                                dnx, dny, dnz,
-                                &newmask1.front(), fbbox1);
-          masklist.push(newmask1);
+          CCTK_INT didit;
           
-          dnx = newbbox2.shape()[0]/newbbox2.stride()[0];
-          dny = newbbox2.shape()[1]/newbbox2.stride()[1];
-          dnz = newbbox2.shape()[2]/newbbox2.stride()[2];
+          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);
           
-          vector<CCTK_INT>  
-            newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0);
+          if (didit == 0) {
           
-          CCTK_FNAME(copy_mask)(nx, ny, nz,
-                                &mask.front(), fbbox,
-                                dnx, dny, dnz,
-                                &newmask2.front(), fbbox2);
-          masklist.push(newmask2);
+            final.push(bb);          
 
-          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());
+            if (verbose) {
+              ostringstream buf;
+              buf << "todo loop. Box pushed to final: " 
+                  << endl << bb;
+              CCTK_INFO(buf.str().c_str());
+            }
           }
-        }
-        else {
-          CCTK_WARN(0, "The fortran routine must be confused.");
-        }
-        
-      } // loop over todo vector (boxes needing to be done).
-      
+          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.");
+          }
+          
+        } // 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;
@@ -450,7 +479,7 @@ namespace CarpetAdaptiveRegrid {
         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);
+          ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel);
         rbbox newbbcoord(lo, hi, str);
         
         if (veryverbose) {
@@ -459,15 +488,15 @@ namespace CarpetAdaptiveRegrid {
               << 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
@@ -475,9 +504,9 @@ namespace CarpetAdaptiveRegrid {
              &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];
@@ -487,30 +516,33 @@ namespace CarpetAdaptiveRegrid {
             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];
-
+          
           newbbcoord = rbbox(lo, up, str);
         }
-        if (veryverbose) {
+        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());
         }
-
+        
         // Convert back to integer coordinates
-
-        ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel);
-        ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel);
+        // 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;
-
+        
         // Check that the width is sufficient
         // This can only be too small if the domain was clipped
         for (int d=0; d < dim; ++d) {
@@ -532,28 +564,27 @@ namespace CarpetAdaptiveRegrid {
             }
           }
         }
-
+        
         ibbox newbb(ilo, ihi, istr);          
-
-        if (veryverbose) {
+        
+        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);
       }
+      
+      
       // FIXME: check if the newbbs is really different from the current bbs
       //        if not, set do_recompose = 0
       bbs = newbbs;
-
-      // make multiprocessor aware
-      gh::cprocs ps;
-      SplitRegions (cctkGH, bbs, obs, ps);
-    
-
+      
+      // Set local bbss
+      
       if (bbss.size() < reflevel+2) {
         if (verbose) {
           CCTK_INFO("Adding new refinement level");
@@ -562,12 +593,19 @@ namespace CarpetAdaptiveRegrid {
         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;
     
-    } // should_regrid?
+    } // did_regrid?
     else
     {
       if (bbss.size() > reflevel+1) {
@@ -577,8 +615,11 @@ namespace CarpetAdaptiveRegrid {
       }
       bbss.resize(reflevel+1);
       obss.resize(reflevel+1);
+      // Set local bbsss
+      MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+      
       pss.resize(reflevel+1);
-
+      
       do_recompose = 1;
     }