Re-organise some refluxing internals

Re-organise some of the internal details of refluxing. Refluxing
requires restricting fluxes from fine to coarse grids. Previously,
Carpet would internally store bounding boxes that were offset by 1/2
grid point, and adjust (correct) these boxes in various places in a
rather ad-hoc manner. This is now cleaned up.

Remove the (unused) "accumulate" operator, and the (unused)
"accumulate" prolongation type.

16 files changed, 102 insertions, 319 deletions

diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index c5279c3f7..3f7618e02 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -2169,8 +2169,6 @@ namespace Carpet {
       return op_restrict;
     } else if (CCTK_Equals(prolong_string, "copy")) {
       return op_copy;
-    } else if (CCTK_Equals(prolong_string, "accumulate")) {
-      return op_accumulate;
     } else if (CCTK_Equals(prolong_string, "Lagrange")) {
       return op_Lagrange;
     } else if (CCTK_Equals(prolong_string, "ENO")) {
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index 552c09130..5503936c8 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -374,8 +374,7 @@ namespace Carpet {
       if (groupdata.AT(group).transport_operator != op_none and
           groupdata.AT(group).transport_operator != op_sync and
           groupdata.AT(group).transport_operator != op_restrict and
-          groupdata.AT(group).transport_operator != op_copy and
-          groupdata.AT(group).transport_operator != op_accumulate)
+          groupdata.AT(group).transport_operator != op_copy)
       {
         if (groupdata.AT(group).activetimelevels.AT(ml).AT(rl) != 0 and
             (groupdata.AT(group).activetimelevels.AT(ml).AT(rl) <
diff --git a/Carpet/CarpetLib/src/accumulate_3d.cc b/Carpet/CarpetLib/src/accumulate_3d.cc
deleted file mode 100644
index 084ad29f2..000000000
--- a/Carpet/CarpetLib/src/accumulate_3d.cc
+++ /dev/null
@@ -1,153 +0,0 @@
-#include <cctk.h>
-#include <cctk_Parameters.h>
-
-#include <algorithm>
-#include <cassert>
-#include <cmath>
-#include <cstdlib>
-#include <iostream>
-
-#include "gdata.hh"
-#include "operator_prototypes_3d.hh"
-#include "typeprops.hh"
-
-using namespace std;
-
-
-
-namespace CarpetLib {
-
-
-  
-#define SRCIND3(i,j,k)                                  \
-  index3 (srcioff + (i), srcjoff + (j), srckoff + (k),  \
-          srciext, srcjext, srckext)
-#define DSTIND3(i,j,k)                                  \
-  index3 (dstioff + (i), dstjoff + (j), dstkoff + (k),  \
-          dstiext, dstjext, dstkext)
-  
-  
-  
-  template <typename T>
-  void
-  accumulate_3d (T const * restrict const src,
-                 ivect3 const & restrict srcext,
-                 T * restrict const dst,
-                 ivect3 const & restrict dstext,
-                 ibbox3 const & restrict srcbbox,
-                 ibbox3 const & restrict dstbbox,
-                 ibbox3 const & restrict,
-                 ibbox3 const & restrict regbbox,
-                 void * extraargs)
-  {
-    assert (not extraargs);
-    
-    if (any (srcbbox.stride() != regbbox.stride() or
-             dstbbox.stride() != regbbox.stride()))
-    {
-#pragma omp critical
-      {
-        cout << "accumulate_3d.cc:" << endl
-             << "srcbbox=" << srcbbox << endl
-             << "dstbbox=" << dstbbox << endl
-             << "regbbox=" << regbbox << endl;
-        CCTK_WARN (0, "Internal error: strides disagree");
-      }
-    }
-    
-    if (any (srcbbox.stride() != dstbbox.stride())) {
-      CCTK_WARN (0, "Internal error: strides disagree");
-    }
-    
-    // This could be handled, but is likely to point to an error
-    // elsewhere
-    if (regbbox.empty()) {
-#pragma omp critical
-      CCTK_WARN (0, "Internal error: region extent is empty");
-    }
-    
-    if (not regbbox.is_contained_in(srcbbox) or
-        not regbbox.is_contained_in(dstbbox))
-    {
-#pragma omp critical
-      {
-        cout << "accumulate_3d.cc:" << endl
-             << "srcbbox=" << srcbbox << endl
-             << "dstbbox=" << dstbbox << endl
-             << "regbbox=" << regbbox << endl;
-        CCTK_WARN (0, "Internal error: region extent is not contained in array extent");
-      }
-    }
-    
-    if (any (srcext != gdata::allocated_memory_shape(srcbbox.shape() / srcbbox.stride()) or
-             dstext != gdata::allocated_memory_shape(dstbbox.shape() / dstbbox.stride())))
-    {
-#pragma omp critical
-      CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes");
-    }
-    
-    
-    
-    ivect3 const regext = regbbox.shape() / regbbox.stride();
-    assert (all ((regbbox.lower() - srcbbox.lower()) % srcbbox.stride() == 0));
-    ivect3 const srcoff = (regbbox.lower() - srcbbox.lower()) / srcbbox.stride();
-    assert (all ((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0));
-    ivect3 const dstoff = (regbbox.lower() - dstbbox.lower()) / dstbbox.stride();
-    
-    
-    
-    ptrdiff_t const srciext = srcext[0];
-    ptrdiff_t const srcjext = srcext[1];
-    ptrdiff_t const srckext = srcext[2];
-    
-    ptrdiff_t const dstiext = dstext[0];
-    ptrdiff_t const dstjext = dstext[1];
-    ptrdiff_t const dstkext = dstext[2];
-    
-    ptrdiff_t const regiext = regext[0];
-    ptrdiff_t const regjext = regext[1];
-    ptrdiff_t const regkext = regext[2];
-    
-    ptrdiff_t const srcioff = srcoff[0];
-    ptrdiff_t const srcjoff = srcoff[1];
-    ptrdiff_t const srckoff = srcoff[2];
-    
-    ptrdiff_t const dstioff = dstoff[0];
-    ptrdiff_t const dstjoff = dstoff[1];
-    ptrdiff_t const dstkoff = dstoff[2];
-    
-    
-    
-    // Loop over region
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
-          
-          dst [DSTIND3(i, j, k)] += src [SRCIND3(i, j, k)];
-          
-        }
-      }
-    }
-    
-  }
-  
-  
-  
-#define TYPECASE(N,T)                                   \
-  template                                              \
-  void                                                  \
-  accumulate_3d (T const * restrict const src,          \
-                 ivect3 const & restrict srcext,        \
-                 T * restrict const dst,                \
-                 ivect3 const & restrict dstext,        \
-                 ibbox3 const & restrict srcbbox,       \
-                 ibbox3 const & restrict dstbbox,       \
-                 ibbox3 const & restrict,               \
-                 ibbox3 const & restrict regbbox,       \
-                 void * extraargs);
-#include "typecase.hh"
-#undef TYPECASE
-  
-  
-  
-} // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 6b28c2e5d..c4f5adfcc 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -539,56 +539,28 @@ copy_from_innerloop (gdata const * const gsrc,
   assert (proc() == src->proc());
   assert (dist::rank() == proc());
   
-  ibbox srcbox, dstbox;
-  switch (cent) {
-  case vertex_centered:
-    srcbox = src->extent();
-    dstbox = this->extent();
-    break;
-  case cell_centered: {
-    ivect const ioff = dstbox.lower() - this->extent().lower();
-    ivect const is_centered = ioff % this->extent().stride() == 0;
-    
-    // Shift bboxes to be face centred if necessary, since all grid
-    // functions are stored as if they were cell-centered
-    // srcbox = src->extent().shift(is_centered-1,2);
-    srcbox = src->extent();
-    dstbox = this->extent().shift(is_centered-1,2);
-    break;
-  }
-  default:
-    assert (0);
-  }
+  ibbox const& srcbox = src->extent();
+  ibbox const& dstbox = this->extent();
   
-  if (transport_operator != op_accumulate) {
 #if CARPET_DIM == 3
-    call_operator<T> (& copy_3d,
-                      static_cast <T const *> (src->storage()),
-                      src->shape(),
-                      static_cast <T *> (this->storage()),
-                      this->shape(),
-                      srcbox, dstbox,
-                      srcregbox, dstregbox, (void*)slabinfo);
+  call_operator<T> (& copy_3d,
+                    static_cast <T const *> (src->storage()),
+                    src->shape(),
+                    static_cast <T *> (this->storage()),
+                    this->shape(),
+                    srcbox, dstbox,
+                    srcregbox, dstregbox, (void*)slabinfo);
 #elif CARPET_DIM == 4
-    call_operator<T> (& copy_4d,
-                      static_cast <T const *> (src->storage()),
-                      src->shape(),
-                      static_cast <T *> (this->storage()),
-                      this->shape(),
-                      srcbox, dstbox,
-                      srcregbox, dstregbox, (void*)slabinfo);
+  call_operator<T> (& copy_4d,
+                    static_cast <T const *> (src->storage()),
+                    src->shape(),
+                    static_cast <T *> (this->storage()),
+                    this->shape(),
+                    srcbox, dstbox,
+                    srcregbox, dstregbox, (void*)slabinfo);
 #else
 #  error "Value for CARPET_DIM not supported"
 #endif
-  } else {
-    call_operator<T> (& accumulate_3d,
-                      static_cast <T const *> (src->storage()),
-                      src->shape(),
-                      static_cast <T *> (this->storage()),
-                      this->shape(),
-                      srcbox, dstbox,
-                      srcregbox, dstregbox, (void*)slabinfo);
-  }
 }
 
 
@@ -699,8 +671,7 @@ transfer_p_r (data const * const src,
   } else if (all (src->extent().stride() < this->extent().stride())) {
     // Restrict
     assert (transport_operator != op_sync);
-    assert (not slabinfo);
-    transfer_restrict (src, dstbox, srcbox, order_space);
+    transfer_restrict (src, dstbox, srcbox, slabinfo, order_space);
   } else {
     assert (0);
   }
@@ -1103,6 +1074,7 @@ data <T>::
 transfer_restrict (data const * const src,
                    ibbox const & dstregbox,
                    ibbox const & srcregbox,
+                   islab const * restrict const slabinfo,
                    int const /*order_space*/)
 {
   static Timer total ("restrict");
@@ -1122,6 +1094,7 @@ transfer_restrict (data const * const src,
     // enum centering { vertex_centered, cell_centered };
     switch (cent) {
     case vertex_centered:
+      assert (not slabinfo);
       call_operator<T> (& restrict_3d_rf2,
                         static_cast <T const *> (src->storage()),
                         src->shape(),
@@ -1133,13 +1106,10 @@ transfer_restrict (data const * const src,
       break;
     case cell_centered: {
       assert (all (dstregbox.stride() == this->extent().stride()));
-      ivect const ioff = dstregbox.lower() - this->extent().lower();
-      ivect const is_centered = ioff % this->extent().stride() == 0;
+      ivect const is_centered = slabinfo ? slabinfo->is_centered : 1;
       
-      // Shift bboxes to be face centred if necessary, since all grid
-      // functions are stored as if they were cell-centered
-      ibbox const srcbox = src->extent().shift(is_centered-1,2);
-      ibbox const dstbox = this->extent().shift(is_centered-1,2);
+      ibbox const& srcbox = src->extent();
+      ibbox const& dstbox = this->extent();
       
       if (all(is_centered == ivect(1,1,1))) {
         call_operator<T> (& restrict_3d_cc_rf2,
@@ -1231,6 +1201,7 @@ data <CCTK_INT>::
 transfer_restrict (data const * const /*src*/,
                    ibbox const & /*dstbox*/,
                    ibbox const & /*srcbox*/,
+                   islab const *restrict const /*slabinfo*/,
                    int const /*order_space*/)
 {
   CCTK_WARN (0, "Data type not supported");
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index ddc645361..f7a79e202 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -106,7 +106,7 @@ private:
                        ibbox const & dstbox,
                        ibbox const & srcbox,
                        islab const * slabinfo);
-  
+
   void
   transfer_from_innerloop (vector <gdata const *> const & gsrcs,
                            vector <CCTK_REAL> const & times,
@@ -150,6 +150,7 @@ private:
   transfer_restrict (data const * const src,
                      ibbox const & dstbox,
                      ibbox const & srcbox,
+                     islab const * restrict const slabinfo,
                      int order_space);
   
   void
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 8dac9b6c4..976a0714f 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -987,14 +987,14 @@ regrid (bool const do_init)
           ibset const tmp3 = tmp2.expand(1,1);
           ibset const tmp4 = all_target - tmp3;
           allrestricted = tmp4;
-          cout << "source=" << source << "\n"
-               << "target=" << target << "\n"
-               << "all=" << all << "\n"
-               << "all_target=" << all_target << "\n"
-               << "tmp1=" << tmp1 << "\n"
-               << "tmp2=" << tmp2 << "\n"
-               << "tmp3=" << tmp3 << "\n"
-               << "allrestricted=" << allrestricted << "\n";
+          // cout << "source=" << source << "\n"
+          //      << "target=" << target << "\n"
+          //      << "all=" << all << "\n"
+          //      << "all_target=" << all_target << "\n"
+          //      << "tmp1=" << tmp1 << "\n"
+          //      << "tmp2=" << tmp2 << "\n"
+          //      << "tmp3=" << tmp3 << "\n"
+          //      << "allrestricted=" << allrestricted << "\n";
           break;
         }
         default:
@@ -1240,8 +1240,10 @@ regrid (bool const do_init)
                       recv.expanded_for (box.exterior.shift(-idir, 2));
                     ASSERT_c (send <= box.exterior.shift(-idir, 2),
                               "Refinement restriction: Send region must be contained in exterior");
-                    
-                    sendrecv_pseudoregion_t const preg (send, c, recv, oc);
+                    ibbox const shifted_recv = recv.shift(idir, 2);
+                    ibbox const shifted_send = send.shift(idir, 2);
+                    sendrecv_pseudoregion_t const preg
+                      (shifted_send, c, shifted_recv, oc);
                     cout << "REF ref_refl ml=" << ml << " rl=" << rl << " olc=" << olc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << "\n"
                          << "   preg=" << preg << "\n";
                     (fast_olevel.*fast_ref_refl_sendrecv).push_back (preg);
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 9c79ed7a4..1a1a18693 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -236,8 +236,7 @@ transfer_from (comm_state & state,
     assert (all(dstbox.lower() >= extent().lower()));
     assert (all(dstbox.upper() <= extent().upper()));
     assert (all(dstbox.stride() == extent().stride()));
-    // This is not satisfied for refluxing
-    // assert (all((dstbox.lower() - extent().lower()) % dstbox.stride() == 0));
+    assert (all((dstbox.lower() - extent().lower()) % dstbox.stride() == 0));
   }
   
   if (is_src) {
@@ -246,10 +245,6 @@ transfer_from (comm_state & state,
     for (int t=0; t<(int)srcs.size(); ++t) {
       assert (srcs.AT(t)->proc() == srcproc);
       assert (srcs.AT(t)->has_storage());
-      if (not slabinfo) {
-        assert (all(srcbox.lower() >= srcs.AT(t)->extent().lower()));
-        assert (all(srcbox.upper() <= srcs.AT(t)->extent().upper()));
-      }
     }
   }
   gdata const * const src = is_src ? srcs.AT(0) : NULL;
@@ -292,16 +287,6 @@ transfer_from (comm_state & state,
       if (is_src) {
         // copy the data into the send buffer
         if (interp_on_src) {
-          ivect ioffset (0);
-          if (src->cent == cell_centered) {
-            assert (all (srcbox.stride() == src->extent().stride()));
-            ivect const ioff = srcbox.lower() - src->extent().lower();
-            ivect const is_centered = ioff % src->extent().stride() == 0;
-            ioffset = 1 - is_centered;
-          }
-          ibbox const bufbox = dstbox.shift(ioffset, 2);
-          assert (prod(allocated_memory_shape(bufbox)) ==
-                  prod(allocated_memory_shape(dstbox)));
           size_t const sendbufsize =
             src->c_datatype_size() * prod(allocated_memory_shape(dstbox));
           void * const sendbuf =
@@ -309,7 +294,7 @@ transfer_from (comm_state & state,
                                prod(allocated_memory_shape(dstbox)));
           gdata * const buf =
             src->make_typed (src->varindex, src->cent, src->transport_operator);
-          buf->allocate (bufbox, srcproc, sendbuf, sendbufsize);
+          buf->allocate (dstbox, srcproc, sendbuf, sendbufsize);
           buf->transfer_from_innerloop
             (srcs, times, dstbox, srcbox, slabinfo,
              time, order_space, order_time);
@@ -327,6 +312,9 @@ transfer_from (comm_state & state,
               src->make_typed (src->varindex, src->cent,
                                src->transport_operator);
             buf->allocate (srcbox, srcproc, sendbuf, sendbufsize);
+            assert (buf->extent().is_aligned_with (this->extent()));
+            assert (srcbox.is_aligned_with (this->extent()));
+            assert (dstbox.is_aligned_with (buf->extent()));
             buf->copy_from_innerloop (srcs.AT(tl), srcbox, srcbox, NULL);
             delete buf;
             state.commit_send_space (src->c_datatype(), dstproc,
@@ -362,13 +350,6 @@ transfer_from (comm_state & state,
           copy_from_innerloop (buf, dstbox, dstbox, NULL);
           delete buf;
         } else {
-          if (cent == cell_centered) {
-            assert (all (dstbox.stride() == this->extent().stride()));
-            ivect const ioff = dstbox.lower() - this->extent().lower();
-            ivect const is_centered = ioff % this->extent().stride() == 0;
-            ivect const ioffset = not is_centered;
-            assert (all (ioffset == 0));
-          }
           gdata const * const null = NULL;
           vector <gdata const *> bufs (ntimelevels, null);
           vector <CCTK_REAL> timebuf (ntimelevels);
@@ -421,7 +402,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
   CCTK_REAL const max_time = * max_element (times.begin(), times.end());
   // TODO: Use a real delta-time from somewhere instead of 1.0
   CCTK_REAL const some_time = abs (min_time) + abs (max_time) + 1.0;
-  if (op != op_copy and op != op_accumulate) {
+  if (op != op_copy) {
     if (time < min_time - eps * some_time or
         time > max_time + eps * some_time)
     {
@@ -449,7 +430,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
     }
   }
   if (timelevel == -1) {
-    if (op == op_copy or op == op_accumulate) {
+    if (op == op_copy) {
       timelevel = 0;
     }
   }
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index 1c903b4f7..1b7e4f203 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -25,8 +25,11 @@ using namespace std;
 // Slabbing description
 template <typename T, int D>
 struct slab {
-  slab() {};
-  vect<T,D> offset;    // dst[ipos] = src[ipos + offset * box.stride];
+  // For periodic boundary conditions:
+  vect<T,D> offset;     // dst[ipos] = src[ipos + offset * box.stride]
+  // For refluxing:
+  vect<T,D> is_centered;
+  slab(): offset(0), is_centered(1) {}
 };
 typedef slab<int,dim> islab;
 
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 9908f0b0b..35bb913c8 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -192,7 +192,8 @@ void ggf::recompose_fill (comm_state & state, int const rl,
       // Initialise from a coarser level of the new hierarchy, where
       // possible
       if (rl > 0) {
-        if (transport_operator != op_none and transport_operator != op_sync and
+        if (transport_operator != op_none and
+            transport_operator != op_sync and
             transport_operator != op_restrict)
         {
           int const numtl = timelevels (ml, rl);
@@ -346,26 +347,6 @@ sync_all (comm_state & state,
 
 
 
-// Accumulate the boundaries of all components
-void
-ggf::
-accumulate_all (comm_state & state,
-                int const tl, int const rl, int const ml)
-{
-  assert (transport_operator == op_accumulate);
-  // Accumulate
-  static Timer timer ("accumulate_all");
-  timer.start ();
-  transfer_from_all (state,
-                     tl,rl,ml,
-                     & dh::fast_dboxes::fast_sync_sendrecv,
-                     tl,rl,ml,
-                     false, true);
-  timer.stop (0);
-}
-
-
-
 // Prolongate the boundaries of all components
 void
 ggf::
@@ -375,13 +356,14 @@ ref_bnd_prolongate_all (comm_state & state,
 {
   // Interpolate
   assert (rl>=1);
-  if (transport_operator == op_none or transport_operator == op_sync or
+  if (transport_operator == op_none or
+      transport_operator == op_sync or
       transport_operator == op_restrict)
     return;
   vector<int> tl2s;
   static Timer timer ("ref_bnd_prolongate_all");
   timer.start ();
-  if (transport_operator != op_copy and transport_operator != op_accumulate) {
+  if (transport_operator != op_copy) {
     // Interpolation in time
     if (not (timelevels(ml,rl) >= prolongation_order_time+1)) {
       char * const fullname = CCTK_FullName (varindex);
@@ -487,7 +469,8 @@ ref_prolongate_all (comm_state & state,
                     CCTK_REAL const time)
 {
   assert (rl>=1);
-  if (transport_operator == op_none or transport_operator == op_sync or
+  if (transport_operator == op_none or
+      transport_operator == op_sync or
       transport_operator == op_restrict)
     return;
   static Timer timer ("ref_prolongate_all");
@@ -514,17 +497,20 @@ ref_reflux_all (comm_state & state,
                 int const tl, int const rl, int const ml,
                 int const dir, int const face)
 {
-  // Ignore the transport operator
   static Timer timer ("ref_reflux_all");
   timer.start ();
   // Require same times
   static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
   assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
           1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
+  islab slabinfo;
+  slabinfo.is_centered = 1 - ivect::dir(dir);
   transfer_from_all (state,
                      tl,rl,ml,
                      dh::fast_dboxes::fast_ref_refl_sendrecv[dir][face],
-                     tl,rl+1,ml);
+                     tl,rl+1,ml,
+                     false, false,
+                     &slabinfo);
   timer.stop (0);
 }
 
@@ -539,7 +525,8 @@ transfer_from_all (comm_state & state,
                    vector<int> const & tl2s, int const rl2, int const ml2,
                    CCTK_REAL const & time,
                    bool const use_old_storage,
-                   bool const flip_send_recv)
+                   bool const flip_send_recv,
+                   islab const *restrict const slabinfo)
 {
   assert (rl1>=0 and rl1<h.reflevels());
   assert (ml1>=0 and ml1<h.mglevels());
@@ -579,11 +566,9 @@ transfer_from_all (comm_state & state,
   // Interpolation orders
   assert (transport_operator != op_none);
   int const pos =
-    transport_operator == op_copy or transport_operator == op_accumulate
-    ? 0 : d.prolongation_orders_space.AT(rl2);
+    transport_operator == op_copy ? 0 : d.prolongation_orders_space.AT(rl2);
   int const pot =
-    transport_operator == op_copy or transport_operator == op_accumulate
-    ? 0 : prolongation_order_time;
+    transport_operator == op_copy ? 0 : prolongation_order_time;
   
   vector<const gdata*> gsrcs(tl2s.size());
   
@@ -626,7 +611,7 @@ transfer_from_all (comm_state & state,
     }
     
     dst->transfer_from
-      (state, gsrcs, times, recv, send, NULL, p1, p2 , time, pos, pot);
+      (state, gsrcs, times, recv, send, slabinfo, p1, p2 , time, pos, pot);
   }
   
   total.stop (0);
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index db2b93f8e..74ebd15ae 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -124,9 +124,6 @@ public:
   // Synchronise the boundaries of a component
   void sync_all (comm_state& state, int tl, int rl, int ml);
 
-  // Accumulate from boundaries of a component
-  void accumulate_all (comm_state& state, int tl, int rl, int ml);
-
   // Prolongate the boundaries of a component
   void ref_bnd_prolongate_all (comm_state& state,
                                int tl, int rl, int ml, CCTK_REAL time);
@@ -172,7 +169,8 @@ protected:
                      vector<int> const & tl2s, int rl2, int ml2,
                      CCTK_REAL const & time,
                      bool use_old_storage = false,
-                     bool flip_send_recv = false);
+                     bool flip_send_recv = false,
+                     islab const *restrict slabinfo = NULL);
   
   void
   transfer_from_all (comm_state & state,
@@ -180,7 +178,8 @@ protected:
                      srpvect const dh::fast_dboxes::* sendrecvs,
                      int tl2, int rl2, int ml2,
                      bool use_old_storage = false,
-                     bool flip_send_recv = false)
+                     bool flip_send_recv = false,
+                     islab const *restrict slabinfo = NULL)
   {
     vector <int> tl2s(1);
     tl2s.AT(0) = tl2;
@@ -191,7 +190,8 @@ protected:
                        tl2s, rl2, ml2,
                        time,
                        use_old_storage,
-                       flip_send_recv);
+                       flip_send_recv,
+                       slabinfo);
   }
 
 
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 45c2b0ee3..4f5c36a47 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -25,7 +25,6 @@ SRCS =	balance.cc				\
 	th.cc					\
 	timestat.cc				\
 	vect.cc					\
-	accumulate_3d.cc			\
 	copy_3d.cc				\
 	copy_4d.cc				\
 	interpolate_3d_2tl.cc			\
diff --git a/Carpet/CarpetLib/src/operator_prototypes_3d.hh b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
index 08d633d18..3cf3cef11 100644
--- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh
+++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
@@ -60,18 +60,6 @@ namespace CarpetLib {
            ibbox3 const & restrict dstregbbox,
            void * extraargs);
   
-  template <typename T>
-  void
-  accumulate_3d (T const * restrict const src,
-                 ivect3 const & restrict srcext,
-                 T * restrict const dst,
-                 ivect3 const & restrict dstext,
-                 ibbox3 const & restrict srcbbox,
-                 ibbox3 const & restrict dstbbox,
-                 ibbox3 const & restrict srcregbbox,
-                 ibbox3 const & restrict dstregbbox,
-                 void * extraargs);
-  
   
   
   template <typename T, int ORDER>
diff --git a/Carpet/CarpetLib/src/operators.hh b/Carpet/CarpetLib/src/operators.hh
index b497430db..af9ae9863 100644
--- a/Carpet/CarpetLib/src/operators.hh
+++ b/Carpet/CarpetLib/src/operators.hh
@@ -12,7 +12,6 @@ enum operator_type
   op_restrict,                  // restrict only, do not prolongate
   op_copy,                      // use simple copying for prolongation
                                 // (needs only one time level)
-  op_accumulate,                // accumulate (sum) into the destination
   op_Lagrange,                  // Lagrange interpolation (standard)
   op_ENO,                       // use ENO stencils (for hydro)
   op_WENO,                      // use WENO stencils (for hydro)
diff --git a/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
index 97047760c..acf27fb95 100644
--- a/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
@@ -135,7 +135,7 @@ namespace CarpetLib {
                       ivect3 const & restrict dstext,
                       ibbox3 const & restrict srcbbox,
                       ibbox3 const & restrict dstbbox,
-                      ibbox3 const & restrict,
+                      ibbox3 const & restrict srcregbbox,
                       ibbox3 const & restrict regbbox,
                       void * extraargs)
   {
@@ -148,9 +148,6 @@ namespace CarpetLib {
     // false: vertex centered, true: cell centered
     ivect const icent (centi, centj, centk);
     assert (all (icent == 0 or icent == 1));
-    bvect const cent (icent);
-    
-    
     
     if (any (srcbbox.stride() >= regbbox.stride() or
              dstbbox.stride() != regbbox.stride()))
@@ -183,22 +180,19 @@ namespace CarpetLib {
     
     
     ivect3 const regext = regbbox.shape() / regbbox.stride();
-    assert (all (srcbbox.stride() % either (cent, 2, 1) == 0));
-    if (not (all ((regbbox.lower() - srcbbox.lower() -
-                   either (cent, srcbbox.stride() / 2, 0)) %
+    assert (all (srcbbox.stride() % 2 == 0));
+    if (not (all ((regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) %
                   srcbbox.stride() == 0)))
     {
       cout << "restrict_3d_vc_rf2.cc\n";
       cout << "regbbox=" << regbbox << "\n";
       cout << "srcbbox=" << srcbbox << "\n";
-      cout << "cent=" << cent << "\n";
+      cout << "icent=" << icent << "\n";
     }
-    assert (all ((regbbox.lower() - srcbbox.lower() -
-                  either (cent, srcbbox.stride() / 2, 0)) %
+    assert (all ((regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) %
                  srcbbox.stride() == 0));
     ivect3 const srcoff =
-      (regbbox.lower() - srcbbox.lower() -
-       either (cent, srcbbox.stride() / 2, 0)) /
+      (regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) /
       srcbbox.stride();
     assert (all ((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0));
     ivect3 const dstoff =
@@ -226,7 +220,9 @@ namespace CarpetLib {
     int const dstjoff = dstoff[1];
     int const dstkoff = dstoff[2];
     
-    size_t const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
+    // size_t const srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0);
+    size_t const srcdi = 1;
+    assert (srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0));
     size_t const srcdj = SRCIND3(0,1,0) - SRCIND3(0,0,0);
     size_t const srcdk = SRCIND3(0,0,1) - SRCIND3(0,0,0);
     
diff --git a/Carpet/CarpetRegrid2/src/property.cc b/Carpet/CarpetRegrid2/src/property.cc
index 404a6e205..c1c8df42c 100644
--- a/Carpet/CarpetRegrid2/src/property.cc
+++ b/Carpet/CarpetRegrid2/src/property.cc
@@ -269,10 +269,10 @@ namespace CarpetRegrid2 {
       // again.
       
       // In detail, we first expand by reffact-1-N points, then expand
-      // to the next coarser grid, then expand back to the current
-      // grid, and expand by N points again. This sequence is correct
-      // for both vertex and cell centred grids, and N is determined
-      // by the number of buffer zones.
+      // (contract???) to the next coarser grid, then expand back to
+      // the current grid, and expand by N points again. This sequence
+      // is correct for both vertex and cell centred grids, and N is
+      // determined by the number of buffer zones.
       
       // N is the number of buffer zones modulo the refinement factor.
       // We cannot shrink the domain (since we cannot shrink
@@ -299,7 +299,17 @@ namespace CarpetRegrid2 {
     if (not snap_to_coarse) return true;
     
     ibset const snapped = snapped_regions (hh, dd, bnd, regions, rl);
-    return regions.AT(rl) == snapped;
+    
+    // We cannot test for equality, since the difference may be
+    // outside of the domain (and hence irrelevant)
+    // return regions.AT(rl) == snapped;
+    
+    // Test whether any part of the difference (i.e. that part of the
+    // level that would be added by snapping) is inside the domain. If
+    // the difference is outside, we can safely ignore it.
+    ibbox const& baseextent = hh.baseextent(0,rl);
+    ibset const difference = snapped - regions.AT(rl);
+    return (difference & baseextent).empty();
   }
   
   void snap_coarse::
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index d656224c2..752d600b7 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -367,6 +367,7 @@ namespace CarpetRegrid2 {
       // Enforce properties on this level
       for (int count=0;; ++count) {
         bool done_enforcing = true;
+        ibset const old_regions = regions.at(rl);
         for (vector<property*>::iterator
                pi = properties.begin(); pi != properties.end(); ++ pi)
         {
@@ -379,8 +380,11 @@ namespace CarpetRegrid2 {
         
         if (done_enforcing) break;
         
+        if (regions.at(rl) == old_regions) {
+          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties (not making any progress); giving up");
+        }
         if (count == 10) {
-          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties; giving up");
+          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties (maximum number of iterations reached); giving up");
         }
         if (count != 0) {
           // This may not be true. However, the previous version of