CarpetLib: Correct prolongation stencils

Correct cell-centered prolongation stencils.
Correct calculation of stencil radii.
Add much more self checking.

2 files changed, 299 insertions, 68 deletions

diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
index 1dd73dfd5..b16f68e62 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
@@ -14,6 +14,33 @@ using namespace std;
 
 
 
+//
+// Grid point locations and their indices:
+//
+// global   0   4      12      20      28   |
+// local    |   0       1       2       3   |
+//          |   C       C       C       C   |
+//          | f   f   f   f   f   f   f   f |
+// local    | 0   1   2   3   4   5   6   7 |
+// global   0 2   6  10  14  18  22  26  30 |
+//
+// Interpolation with even interpolation order:
+//    offset zero (fine index = 2 * coarse index):
+//       [1]
+//    offset one (fine index = 2 * coarse index + 1):
+//       [1]
+// Interpolation with odd interpolation order:
+//    offset zero:
+//       [1 1 0]/2
+//    offset one:
+//       [0 1 1]/2
+// The centres of these stencils are located at the coarse grid point
+// corresponding to the fine grid point (the interpolation target)
+// minus the offset. Example: fine grid 8 -> coarse grid 8, fine grid
+// 12 -> also coarse grid 8.
+
+
+
 namespace CarpetLib {
   
 
@@ -25,28 +52,37 @@ namespace CarpetLib {
   
     
   namespace coeffs_3d_cc_rf2 {
-
+    
   // 1D interpolation coefficients
   
   template<typename RT, int ORDER>
-  struct coeffs1d {
+  struct coeffs1dc {
     static const RT coeffs[];
+  };
+  
+  template<typename RT, int ORDER, int di>
+  struct coeffs1d {
+    typedef coeffs1dc<RT,ORDER> coeffs_t;
     
     static ptrdiff_t const ncoeffs = ORDER+1;
-    static ptrdiff_t const imin    = - ncoeffs/2 + 1;
+    static ptrdiff_t const imin    = - ncoeffs/2 + (ORDER % 2 != 0 ? di : 0);
     static ptrdiff_t const imax    = imin + ncoeffs;
     
     static inline
     RT
-    get (int const di, int const i)
+    get (ptrdiff_t const i)
     {
-      ptrdiff_t const j = di == 0 ? i - imin : imax-1 - (i-imin);
+      static_assert (ncoeffs == sizeof coeffs_t::coeffs / sizeof *coeffs_t::coeffs,
+                     "coefficient array has wrong size");
+      static_assert (di==0 or di==1, "di must be 0 or 1");
+#ifdef CARPET_DEBUG
+      assert (i>=imin and i<imax);
+#endif
+      ptrdiff_t const j = di==0 ? i-imin : ncoeffs-1 - (i-imin);
 #ifdef CARPET_DEBUG
-      assert (di == 0 or di == 1);
-      assert (ncoeffs == sizeof coeffs / sizeof *coeffs);
       assert (j>=0 and j<ncoeffs);
 #endif
-      return coeffs[j];
+      return coeffs_t::coeffs[j];
     }
     
     // Test coefficients
@@ -58,35 +94,33 @@ namespace CarpetLib {
       if (tested) return;
       tested = true;
       
-      assert (ncoeffs == sizeof coeffs / sizeof *coeffs);
+      static_assert (ncoeffs == sizeof coeffs_t::coeffs / sizeof *coeffs_t::coeffs,
+                     "coefficient array has wrong size");
       
-      // Test all orders and offsets
+      // Test all orders
       bool error = false;
       for (int n=0; n<=ORDER; ++n) {
-        for (int di=0; di<2; ++di) {
-          RT res = RT(0.0);
-          for (int i=imin; i<imax; ++i) {
-            CCTK_REAL const dx = ORDER % 2;
-            RT const x = RT(CCTK_REAL(i) - (di==0 ? 0.75 : 1.25 - dx));
-            RT const y = ipow (x, n);
-            res += get(di,i) * y;
-          }
-          RT const x0 = RT(0.0);
-          RT const y0 = ipow (x0, n);
-          if (not (good::abs (res - y0) < 1.0e-12)) {
-            RT rt;
-            ostringstream buf;
-            buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
-                << "   RT=" << typestring(rt) << "\n"
-                << "   ORDER=" << ORDER << "\n"
-                << "   n=" << n << "\n"
-                << "   di=" << di << "\n"
-                << "   y0=" << y0 << ", res=" << res;
-            CCTK_WARN (CCTK_WARN_ALERT, buf.str().c_str());
-            error = true;
-          }
-        } // for di
-      }   // for n
+        RT res = RT(0.0);
+        for (ptrdiff_t i=imin; i<imax; ++i) {
+          RT const x = RT(CCTK_REAL(i) + 0.5);
+          RT const y = ipow (x, n);
+          res += get(i) * y;
+        }
+        RT const x0 = RT(0.25) + di * RT(0.5);
+        RT const y0 = ipow (x0, n);
+        if (not (good::abs (res - y0) < 1.0e-12)) {
+          RT rt;
+          ostringstream buf;
+          buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
+              << "   RT=" << typestring(rt) << "\n"
+              << "   ORDER=" << ORDER << "\n"
+              << "   di=" << di << "\n"
+              << "   n=" << n << "\n"
+              << "   y0=" << y0 << ", res=" << res;
+          CCTK_WARN (CCTK_WARN_ALERT, buf.str().c_str());
+          error = true;
+        }
+      } // for n
       if (error) {
         CCTK_WARN (CCTK_WARN_ABORT, "Aborting.");
       }
@@ -94,29 +128,29 @@ namespace CarpetLib {
   };
   
   
-   
+  
 #define TYPECASE(N,RT)                          \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,0>::coeffs[] = {         \
+  const RT coeffs1dc<RT,0>::coeffs[] = {        \
     +1 / RT(1.0)                                \
   };                                            \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,1>::coeffs[] = {         \
+  const RT coeffs1dc<RT,1>::coeffs[] = {        \
     +1 / RT(4.0),                               \
     +3 / RT(4.0)                                \
   };                                            \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,2>::coeffs[] = {         \
+  const RT coeffs1dc<RT,2>::coeffs[] = {        \
     + 5 / RT(32.0),                             \
     +30 / RT(32.0),                             \
     - 3 / RT(32.0)                              \
   };                                            \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,3>::coeffs[] = {         \
+  const RT coeffs1dc<RT,3>::coeffs[] = {        \
     -  5 / RT(128.0),                           \
     + 35 / RT(128.0),                           \
     +105 / RT(128.0),                           \
@@ -124,24 +158,24 @@ namespace CarpetLib {
   };                                            \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,4>::coeffs[] = {         \
+  const RT coeffs1dc<RT,4>::coeffs[] = {        \
     -  45 / RT(2048.0),                         \
     + 420 / RT(2048.0),                         \
     +1890 / RT(2048.0),                         \
-    - 252 / RT(2048.0)                          \
+    - 252 / RT(2048.0),                         \
     +  35 / RT(2048.0)                          \
   };                                            \
                                                 \
   template<>                                    \
-  const RT coeffs1d<RT,5>::coeffs[] = {         \
+  const RT coeffs1dc<RT,5>::coeffs[] = {        \
     +  63 / RT(8192.0),                         \
     - 495 / RT(8192.0),                         \
     +2310 / RT(8192.0),                         \
     +6930 / RT(8192.0),                         \
-    - 639 / RT(8192.0),                         \
+    - 693 / RT(8192.0),                         \
     +  77 / RT(8192.0)                          \
   };
-
+  
 #define CARPET_NO_COMPLEX
 #include "typecase.hh"
 #undef CARPET_NO_COMPLEX
@@ -169,11 +203,12 @@ namespace CarpetLib {
   interp1 (T const * restrict const p,
            size_t const d1)
   {
+    static_assert (di==0 or di==1, "di must be 0 or 1");
     typedef typename typeprops<T>::real RT;
-    typedef coeffs1d<RT,ORDER> coeffs;
+    typedef coeffs1d<RT,ORDER,di> coeffs;
     T res = typeprops<T>::fromreal (0);
     for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) {
-      res += coeffs::get(di,i) * interp0<T,ORDER> (p + i*d1);
+      res += coeffs::get(i) * interp0<T,ORDER> (p + i*d1);
     }
     return res;
   }
@@ -186,11 +221,12 @@ namespace CarpetLib {
            size_t const d1,
            size_t const d2)
   {
+    static_assert (dj==0 or dj==1, "dj must be 0 or 1");
     typedef typename typeprops<T>::real RT;
-    typedef coeffs1d<RT,ORDER> coeffs;
+    typedef coeffs1d<RT,ORDER,dj> coeffs;
     T res = typeprops<T>::fromreal (0);
     for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) {
-      res += coeffs::get(dj,i) * interp1<T,ORDER,di> (p + i*d2, d1);
+      res += coeffs::get(i) * interp1<T,ORDER,di> (p + i*d2, d1);
     }
     return res;
   }
@@ -204,17 +240,79 @@ namespace CarpetLib {
            size_t const d2,
            size_t const d3)
   {
+    static_assert (dk==0 or dk==1, "dk must be 0 or 1");
     typedef typename typeprops<T>::real RT;
-    typedef coeffs1d<RT,ORDER> coeffs;
+    typedef coeffs1d<RT,ORDER,dk> coeffs;
     T res = typeprops<T>::fromreal (0);
     for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) {
-      res += coeffs::get(dk,i) * interp2<T,ORDER,di,dj> (p + i*d3, d1, d2);
+      res += coeffs::get(i) * interp2<T,ORDER,di,dj> (p + i*d3, d1, d2);
     }
     return res;
   }
   
   
   
+  // Check interpolation index ranges
+  template <typename T, int ORDER>
+  static inline
+  void
+  check_indices0 ()
+  {
+  }
+  
+  template <typename T, int ORDER, int di>
+  static inline
+  void
+  check_indices1 (ptrdiff_t const is,
+                  ptrdiff_t const srciext)
+  {
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER,di> coeffs;
+    assert (is + coeffs::imin >= 0);
+    assert (is + coeffs::imax <= srciext);
+    check_indices0<T,ORDER> ();
+#endif
+  }
+  
+  template <typename T, int ORDER, int di, int dj>
+  static inline
+  void
+  check_indices2 (ptrdiff_t const is,
+                  ptrdiff_t const js,
+                  ptrdiff_t const srciext,
+                  ptrdiff_t const srcjext)
+  {
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER,dj> coeffs;
+    assert (js + coeffs::imin >= 0);
+    assert (js + coeffs::imax <= srcjext);
+    check_indices1<T,ORDER,di> (is, srciext);
+#endif
+  }
+  
+  template <typename T, int ORDER, int di, int dj, int dk>
+  static inline
+  void
+  check_indices3 (ptrdiff_t const is,
+                  ptrdiff_t const js,
+                  ptrdiff_t const ks,
+                  ptrdiff_t const srciext,
+                  ptrdiff_t const srcjext,
+                  ptrdiff_t const srckext)
+  {
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER,dk> coeffs;
+    assert (ks + coeffs::imin >= 0);
+    assert (ks + coeffs::imax <= srckext);
+    check_indices2<T,ORDER,di,dj> (is,js, srciext,srcjext);
+#endif
+  }
+  
+  
+  
   template <typename T, int ORDER>
   void
   prolongate_3d_cc_rf2 (T const * restrict const src,
@@ -228,7 +326,8 @@ namespace CarpetLib {
     static_assert (ORDER>=0, "ORDER must be non-negative");
     
     typedef typename typeprops<T>::real RT;
-    coeffs1d<RT,ORDER>::test ();
+    coeffs1d<RT,ORDER,0>::test ();
+    coeffs1d<RT,ORDER,1>::test ();
     
     
     
@@ -262,14 +361,27 @@ namespace CarpetLib {
     
     
     
-    // The name "needoffset" does not make sense for cell centring
-    bvect3 const needoffsetlo = srcoff % reffact2 != 0;
-    bvect3 const needoffsethi = (srcoff + regext - 1) % reffact2 != 0;
-    // This is probably wrong for odd orders
-    ivect3 const offsetlo =
-      ORDER%2!=0 ? ORDER/2 : either (needoffsetlo, ORDER/2+1, ORDER/2);
-    ivect3 const offsethi =
-      ORDER%2!=0 ? ORDER/2 : either (needoffsethi, ORDER/2+1, ORDER/2);
+    // Determine the stencil radius (see diagram at the top of this
+    // file)
+    ivect3 offsetlo, offsethi;
+    {
+      assert (all ((regbbox.upper() - srcbbox.lower() + regbbox.stride() / 2) % regbbox.stride() == 0));
+      ivect3 const srcend = (regbbox.upper() - srcbbox.lower() + regbbox.stride() / 2) / regbbox.stride();
+      ivect const needoffsetlo = srcoff % 2;
+      ivect const needoffsethi = srcend % 2;
+      for (int d=0; d<3; ++d) {
+        if (not needoffsetlo[d]) {
+          offsetlo[d] = - coeffs1d<RT,ORDER,0>::imin;
+        } else {
+          offsetlo[d] = - coeffs1d<RT,ORDER,1>::imin;
+        }
+        if (not needoffsethi[d]) {
+          offsethi[d] = coeffs1d<RT,ORDER,0>::imax - 1;
+        } else {
+          offsethi[d] = coeffs1d<RT,ORDER,1>::imax - 1;
+        }
+      }
+    }
     
     
     
@@ -327,8 +439,8 @@ namespace CarpetLib {
     
     
     // size_t const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
-    assert (SRCIND3(1,0,0) - SRCIND3(0,0,0) == 1);
     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);
     
@@ -366,6 +478,7 @@ namespace CarpetLib {
     
     // kernel
   l8000:
+    check_indices3<T,ORDER,0,0,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,0,0> (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
     i = i+1;
@@ -375,6 +488,7 @@ namespace CarpetLib {
     
     // kernel
   l8001:
+    check_indices3<T,ORDER,1,0,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,0,0> (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
     i = i+1;
@@ -400,6 +514,7 @@ namespace CarpetLib {
     
     // kernel
   l8010:
+    check_indices3<T,ORDER,0,1,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,1,0> (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
     i = i+1;
@@ -409,6 +524,7 @@ namespace CarpetLib {
     
     // kernel
   l8011:
+    check_indices3<T,ORDER,1,1,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,1,0> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -450,6 +566,7 @@ namespace CarpetLib {
     
     // kernel
   l8100:
+    check_indices3<T,ORDER,0,0,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,0,1> (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
     i = i+1;
@@ -459,6 +576,7 @@ namespace CarpetLib {
     
     // kernel
   l8101:
+    check_indices3<T,ORDER,1,0,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,0,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -484,6 +602,7 @@ namespace CarpetLib {
     
     // kernel
   l8110:
+    check_indices3<T,ORDER,0,1,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,1,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -493,6 +612,7 @@ namespace CarpetLib {
     
     // kernel
   l8111:
+    check_indices3<T,ORDER,1,1,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,1,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
diff --git a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
index 4f669d6ea..fb85a9503 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
@@ -43,9 +43,13 @@ namespace CarpetLib {
     RT const&
     get (ptrdiff_t const i)
     {
+      static_assert (ncoeffs == sizeof coeffs / sizeof *coeffs,
+                     "coefficient array has wrong size");
+#ifdef CARPET_DEBUG
+      assert (i>=imin and i<imax);
+#endif
       ptrdiff_t const j = i - imin;
 #ifdef CARPET_DEBUG
-      assert (ncoeffs == sizeof coeffs / sizeof *coeffs);
       assert (j>=0 and j<ncoeffs);
 #endif
       return coeffs[j];
@@ -60,7 +64,8 @@ namespace CarpetLib {
       if (tested) return;
       tested = true;
       
-      assert (ncoeffs == sizeof coeffs / sizeof *coeffs);
+      static_assert (ncoeffs == sizeof coeffs / sizeof *coeffs,
+                     "coefficient array has wrong size");
       
       // Do not test integer operators (they should be disabled
       // anyway)
@@ -68,7 +73,7 @@ namespace CarpetLib {
       
       // Test all orders
       bool error = false;
-      for (ptrdiff_t n=0; n<=ORDER; ++n) {
+      for (int n=0; n<=ORDER; ++n) {
         RT res = RT(0.0);
         for (ptrdiff_t i=imin; i<imax; ++i) {
           RT const x = RT(CCTK_REAL(i) - 0.5);
@@ -363,6 +368,85 @@ namespace CarpetLib {
   
   
   
+  // Check interpolation index ranges
+  template <typename T, int ORDER>
+  static inline
+  void
+  check_indices0 ()
+  {
+  }
+  
+  template <typename T, int ORDER, int di>
+  static inline
+  void
+  check_indices1 (ptrdiff_t const is,
+                  ptrdiff_t const srciext)
+  {
+    static_assert (di==0 or di==1, "di must be 0 or 1");
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER> coeffs;
+    if (di == 0) {
+      assert (is >= 0);
+      assert (is < srciext);
+    } else {
+      assert (is + coeffs::imin >= 0);
+      assert (is + coeffs::imax <= srciext);
+    }
+    check_indices0<T,ORDER> ();
+#endif
+  }
+  
+  template <typename T, int ORDER, int di, int dj>
+  static inline
+  void
+  check_indices2 (ptrdiff_t const is,
+                  ptrdiff_t const js,
+                  ptrdiff_t const srciext,
+                  ptrdiff_t const srcjext)
+  {
+    static_assert (dj==0 or dj==1, "dj must be 0 or 1");
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER> coeffs;
+    if (dj == 0) {
+      assert (js >= 0);
+      assert (js < srcjext);
+    } else {
+      assert (js + coeffs::imin >= 0);
+      assert (js + coeffs::imax <= srcjext);
+    }
+    check_indices1<T,ORDER,di> (is, srciext);
+#endif
+  }
+  
+  template <typename T, int ORDER, int di, int dj, int dk>
+  static inline
+  void
+  check_indices3 (ptrdiff_t const is,
+                  ptrdiff_t const js,
+                  ptrdiff_t const ks,
+                  ptrdiff_t const srciext,
+                  ptrdiff_t const srcjext,
+                  ptrdiff_t const srckext)
+  {
+    static_assert (dk==0 or dk==1, "dk must be 0 or 1");
+#ifdef CARPET_DEBUG
+    typedef typename typeprops<T>::real RT;
+    typedef coeffs1d<RT,ORDER> coeffs;
+    if (dk == 0) {
+      assert (ks >= 0);
+      assert (ks < srckext);
+    } else {
+      assert (ks + coeffs::imin >= 0);
+      assert (ks + coeffs::imax <= srckext);
+    }
+    check_indices2<T,ORDER,di,dj> (is,js, srciext,srcjext);
+#endif
+  }
+  
+  
+  
   template <typename T, int ORDER>
   void
   prolongate_3d_rf2 (T const * restrict const src,
@@ -391,6 +475,16 @@ namespace CarpetLib {
       CCTK_WARN (0, "Internal error: source strides are not twice the destination strides");
     }
     
+    if (any (srcbbox.lower() % srcbbox.stride() != 0)) {
+      CCTK_WARN (0, "Internal error: source bbox is not aligned with vertices");
+    }
+    if (any (dstbbox.lower() % dstbbox.stride() != 0)) {
+      CCTK_WARN (0, "Internal error: destination bbox is not aligned with vertices");
+    }
+    if (any (regbbox.lower() % regbbox.stride() != 0)) {
+      CCTK_WARN (0, "Internal error: prolongation region bbox is not aligned with vertices");
+    }
+    
     // This could be handled, but is likely to point to an error
     // elsewhere
     if (regbbox.empty()) {
@@ -407,16 +501,25 @@ namespace CarpetLib {
     
     
     
-    bvect3 const needoffsetlo = srcoff % reffact2 != 0 or regext > 1;
-    bvect3 const needoffsethi = (srcoff + regext - 1) % reffact2 != 0 or regext > 1;
-    ivect3 const offsetlo = either (needoffsetlo, ORDER/2+1, 0);
-    ivect3 const offsethi = either (needoffsethi, ORDER/2+1, 0);
+    bvect3 const needoffsetlo = srcoff % reffact2 != 0;
+    bvect3 const needoffsethi = (srcoff + regext - 1) % reffact2 != 0;
+    ivect3 const offsetlo =
+      either (needoffsetlo, ORDER/2+1, either (regext > 1, ORDER/2, 0));
+    ivect3 const offsethi =
+      either (needoffsethi, ORDER/2+1, either (regext > 1, ORDER/2, 0));
     
     
     
     if (not regbbox.expand(offsetlo, offsethi).is_contained_in(srcbbox) or
         not regbbox                           .is_contained_in(dstbbox))
     {
+      cerr << "ORDER=" << ORDER << "\n"
+           << "offsetlo=" << offsetlo << "\n"
+           << "offsethi=" << offsethi << "\n"
+           << "regbbox=" << regbbox << "\n"
+           << "dstbbox=" << dstbbox << "\n"
+           << "regbbox.expand=" << regbbox.expand(offsetlo, offsethi) << "\n"
+           << "srcbbox=" << srcbbox << "\n";
       CCTK_WARN (0, "Internal error: region extent is not contained in array extent");
     }
     
@@ -461,8 +564,8 @@ namespace CarpetLib {
     
     
     // size_t const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
-    assert (SRCIND3(1,0,0) - SRCIND3(0,0,0) == 1);
     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);
     
@@ -500,6 +603,7 @@ namespace CarpetLib {
     
     // kernel
   l8000:
+    check_indices3<T,ORDER,0,0,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,0,0> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -509,6 +613,7 @@ namespace CarpetLib {
     
     // kernel
   l8001:
+    check_indices3<T,ORDER,1,0,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,0,0> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -534,6 +639,7 @@ namespace CarpetLib {
     
     // kernel
   l8010:
+    check_indices3<T,ORDER,0,1,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,1,0> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -543,6 +649,7 @@ namespace CarpetLib {
     
     // kernel
   l8011:
+    check_indices3<T,ORDER,1,1,0> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,1,0> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -584,6 +691,7 @@ namespace CarpetLib {
     
     // kernel
   l8100:
+    check_indices3<T,ORDER,0,0,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,0,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -593,6 +701,7 @@ namespace CarpetLib {
     
     // kernel
   l8101:
+    check_indices3<T,ORDER,1,0,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,0,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -618,6 +727,7 @@ namespace CarpetLib {
     
     // kernel
   l8110:
+    check_indices3<T,ORDER,0,1,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,0,1,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;
@@ -627,6 +737,7 @@ namespace CarpetLib {
     
     // kernel
   l8111:
+    check_indices3<T,ORDER,1,1,1> (is,js,ks, srciext, srcjext, srckext);
     dst[DSTIND3(id,jd,kd)] =
       interp3<T,ORDER,1,1,1> (& src[SRCIND3(is,js,ks)], srcdi,srcdj,srcdk);
     i = i+1;