CarpetLib: clean up higher order restriciton some

* rename controlling parameter to use_higher_order_restriction
* introduce parameter restriction_order_space to control which operator
  is used (currently order 1 and 3 are suppoted)
* include some comments on what the operator does
* change the way the restrictable region is computed in dh.cc/regrid to
  be based on exterior.shrink(stencil_width) rather that the interior

4 files changed, 53 insertions, 16 deletions

diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index b0a825171..a259fff2d 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -33,10 +33,16 @@ BOOLEAN use_dgfe "Use DGFE operators instead of Lagrange operators"
 
 
 
-BOOLEAN use_cc_o3 "Use third order cell centered restriction operators instead of first order"
+BOOLEAN use_higher_order_restriction "Use third order cell centered restriction operators instead of first order"
 {
 } "no"
 
+INT restriction_order_space "Order of restriction operator to use with use_higher_order_restriction"
+{
+  1 :: "linear interpolation, this is Carpet's original implementation"
+  3 :: "third order accurate restriction for grid functions where prolongation is not (W)ENO"
+} 3
+
 
 
 BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=always
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 9a26eec8a..ba1a438ce 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1141,17 +1141,35 @@ transfer_restrict (data const * const src,
                                     srcregbox, dstregbox, NULL);
           break;
         }
-        if (use_cc_o3 and
+        if (use_higher_order_restriction and
             transport_operator != op_WENO and 
             transport_operator != op_ENO) { // HACK
-          // Don't use call_operator, because we parallelise ourselves
-          restrict_3d_cc_o3_rf2(static_cast <T const *> (src->storage()),
-                             src->shape(),
-                             static_cast <T *> (this->storage()),
-                             this->shape(),
-                             srcbox,
-                             dstbox,
-                             srcregbox, dstregbox, NULL);
+	  switch (restriction_order_space) {
+	    case 1:
+            // Don't use call_operator, because we parallelise ourselves
+            restrict_3d_cc_rf2(static_cast <T const *> (src->storage()),
+                               src->shape(),
+                               static_cast <T *> (this->storage()),
+                               this->shape(),
+                               srcbox,
+                               dstbox,
+                               srcregbox, dstregbox, NULL);
+	    break;
+	    case 3:
+            // Don't use call_operator, because we parallelise ourselves
+            restrict_3d_cc_o3_rf2(static_cast <T const *> (src->storage()),
+                               src->shape(),
+                               static_cast <T *> (this->storage()),
+                               this->shape(),
+                               srcbox,
+                               dstbox,
+                               srcregbox, dstregbox, NULL);
+	    break;
+	    default:
+	    CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+			"There is no restriction stencil with restriction_order_space==%d", restriction_order_space);
+	    break;
+	  }
           break;
         }
         // Don't use call_operator, because we parallelise ourselves
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 572fe2e75..1380a762a 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1034,7 +1034,7 @@ regrid (bool const do_init)
           full_dboxes const& obox = full_olevel.AT(oc);
           
           ibset needrecv = allrestricted & obox.owned;
-          if(use_cc_o3) {
+          if(use_higher_order_restriction) {
             // NOTE: change in behaviour (affects only outer boundaries I think)!!!!
             // NOTE: b/c of this we need a low-level sync after the restrict
             needrecv = allrestricted & obox.interior;
@@ -1045,10 +1045,14 @@ regrid (bool const do_init)
           for (int c = 0; c < h.components(rl); ++ c) {
             full_dboxes const& box = full_level.AT(c);
             
-            // HORRIBLE HACK
-            ibbox const contracted_exterior =
-              use_cc_o3 ? box.interior.expand(ivect(int(h.refcent==cell_centered))).contracted_for(odomext) :
-              box.exterior.contracted_for(odomext);
+            // If we make a mistake expanding the domain of dependence here, it
+            // _should_ be caught be by the expand()ed is_contained_in(srcbbox)
+            // test in the actual operator. 
+            int const shrink_by =
+              use_higher_order_restriction and (h.refcent == cell_centered) ?
+              restriction_order_space/2  : 0;
+            ibbox const contracted_exterior = 
+              box.exterior.expand(ivect(-shrink_by)).contracted_for(odomext);
             ibset const ovlp = needrecv & contracted_exterior;
             
             for (ibset::const_iterator
@@ -1058,7 +1062,7 @@ regrid (bool const do_init)
               ibbox const send = recv.expanded_for(box.exterior);
               ASSERT_c (send <= box.exterior,
                         "Refinement restriction: Send region must be contained in exterior");
-              if(use_cc_o3) {
+              if(use_higher_order_restriction) {
                 ASSERT_c (send <= box.interior.expand(ivect(int(h.refcent==cell_centered))),
                           "Refinement restriction: Send region must be contained in interior");
               }
diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
index c8ecc040f..991734c85 100644
--- a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
@@ -25,6 +25,15 @@ namespace CarpetLib {
           dstiext, dstjext, dstkext)
   
   
+// This operator offers third-order accurate restriction operators for cell
+// centered grids when use_higher_order_restriction is set. This interpolation
+// is done for samples at the cell centres (so this is not a ppm scheme or
+// anything like that).  It really only fits a polynomial of the form
+//    f(x,y,z) = \sum_{i,j,k=0}^3 a_{i,j,k} x^i y^j z^k
+// to the fine cells and evaluates at x=y=z=0. So it is good for the metric,
+// but bad for matter (since it will destroy the conservation).  Because of
+// this it should not be used for grid functions whose transport operator is
+// not WENO or ENO which hopefully excludes all matter variables.
   
   template <typename T>
   void