CarpetLib: Correct calculation of required stencil sizes

Correct the calculation of required stencil sizes in the cell-centred prolongation operators.

1 files changed, 15 insertions, 4 deletions

diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
index 3bd9e6aad..1dd73dfd5 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
@@ -6,6 +6,7 @@
 #include <cmath>
 #include <cstdlib>
 
+#include "gdata.hh"
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
@@ -261,22 +262,32 @@ 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 = either (needoffsetlo, ORDER/2+1, ORDER/2);
-    ivect3 const offsethi = either (needoffsethi, ORDER/2+1, ORDER/2);
+    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);
     
     
     
     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");
     }
     
-    if (any (srcext != srcbbox.shape() / srcbbox.stride() or
-             dstext != dstbbox.shape() / dstbbox.stride()))
+    if (any (srcext != gdata::allocated_memory_shape(srcbbox.shape() / srcbbox.stride()) or
+             dstext != gdata::allocated_memory_shape(dstbbox.shape() / dstbbox.stride())))
     {
       CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes");
     }