Update for cartoon.

7 files changed, 86 insertions, 28 deletions

diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index f88abdb51..1f007f4d1 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -14,6 +14,7 @@
 #include <string>
 #include <vector>
 
+#include <unistd.h>
 #include <cctk.h>
 #include <cctk_Parameters.h>
 #include <util_ErrorCodes.h>
@@ -199,6 +200,10 @@ namespace Carpet {
            cGH * const cctkGH)
   {
     DECLARE_CCTK_PARAMETERS;
+
+    int i = 1;
+    //while (i)
+    //    usleep(1000000);
     
     // Some statistics
     {
@@ -2780,9 +2785,9 @@ namespace Carpet {
         vdd.AT(m)->prolongation_orders_space.AT(rl);
       int const prolongation_stencil_size =
         vdd.AT(m)->prolongation_stencil_size(rl);
-      int const min_nghosts =
-        ((prolongation_stencil_size + refinement_factor - 1)
-         / (refinement_factor - 1));
+      int const min_nghosts = 0;
+        //((prolongation_stencil_size + refinement_factor - 1)
+        // / (refinement_factor - 1));
       int const min_nghosts_restrict =
         restriction_order_space / 2;
       if (any (any (ghosts.AT(rl) < i2vect (min_nghosts)))) {
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 3418a932c..0e1474e21 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -367,6 +367,7 @@ namespace CarpetInterp {
     int const minrl = want_global_mode ? 0 : reflevel;
     int const maxrl = want_global_mode ?
       arrdata.AT(coord_group).AT(0).hh->reflevels() : reflevel+1;
+    //int maxrl = reflevel + 1;
 
     // Find maximum number of components over all levels and maps
     int maxncomps = 0;
@@ -1339,7 +1340,9 @@ namespace CarpetInterp {
       
       // Number of neccessary time levels
       // CCTK_REAL const level_time = cctkGH->cctk_time;
-      CCTK_REAL const level_time = tt->get_time(mglevel, rl, tl);
+      //CCTK_REAL const level_time = tt->get_time(mglevel, rl, tl);
+      double get_time = tt->get_time(mglevel, rl, tl);
+      CCTK_REAL const level_time = (rl == reflevel && tl == 0) ? cctkGH->cctk_time : tt->get_time(mglevel, rl, tl);
       vector<int> num_tl (N_input_arrays, 0);
       vector<bool> need_time_interp (N_output_arrays);
       for (int m=0; m<N_output_arrays; ++m) {
@@ -1405,11 +1408,11 @@ namespace CarpetInterp {
   class InterpolationTimes : private vector<CCTK_REAL>
   {
   public:
-    InterpolationTimes (int const rl, int const num_timelevels_ )
+    InterpolationTimes (const cGH *cctkGH, int reflevel, int const rl, int const num_timelevels_ )
       : vector<CCTK_REAL> (num_timelevels_ )
     {
       for (int tl=0; tl<num_timelevels_; ++tl) {
-        at(tl) = tt->get_time (mglevel, rl, tl);
+        at(tl) = (rl == reflevel && tl == 0) ? cctkGH->cctk_time : tt->get_time(mglevel, rl, tl);
       }
     }
 
@@ -1767,7 +1770,7 @@ namespace CarpetInterp {
         int const interp_num_tl =
           interp_num_time_levels.AT(n) > 0 ?
           interp_num_time_levels.AT(n) : num_tl.AT(n);
-        const InterpolationTimes times (rl, interp_num_tl);
+        const InterpolationTimes times (cctkGH, reflevel, rl, interp_num_tl);
         const InterpolationWeights tfacs (deriv_order, times, current_time,
                                           delta_time);
         
diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh
index 5dea7f9d3..b2eb6de1f 100644
--- a/Carpet/CarpetLib/src/bbox.hh
+++ b/Carpet/CarpetLib/src/bbox.hh
@@ -43,13 +43,13 @@ class bbox {
   // Fields
   
   /** Bounding box bounds and stride.  The bounds are inclusive.  */
-  vect<T,D> _lower, _upper, _stride;
   
   // Consistency checks
   
   void assert_bbox_limits () const CCTK_MEMBER_ATTRIBUTE_PURE;
   
 public:
+  vect<T,D> _lower, _upper, _stride;
   
   // Constructors
   
diff --git a/Carpet/CarpetLib/src/bboxset1.hh b/Carpet/CarpetLib/src/bboxset1.hh
index db760e5ae..e268db588 100644
--- a/Carpet/CarpetLib/src/bboxset1.hh
+++ b/Carpet/CarpetLib/src/bboxset1.hh
@@ -85,7 +85,6 @@ class bboxset {
 #endif
   
   // Fields
-  bset bs;
   // Invariant:
   // All bboxes have the same stride.
   // No bbox is empty.
@@ -94,6 +93,7 @@ class bboxset {
   bool skip_normalize;
   
 public:
+  bset bs;
   
   // Constructors
   bboxset ();                   // cost: O(1)
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index b82337685..7f914e1f6 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,5 +1,6 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
+#include <unistd.h>
 
 #include <cassert>
 #include <cstddef>
@@ -772,7 +773,10 @@ regrid (bool const do_init)
           
           ibset oneedrecv = obox.active;
           
-          i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
+          int ssize = prolongation_stencil_size(rl);
+          ivect ss = ivect(ssize, 0, ssize);
+          i2vect const stencil_size = i2vect (ss, ss);
+          //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
           
           ibset const expanded_active (box.active.expanded_for (obox.interior));
           ibset const ovlp = oneedrecv & expanded_active;
@@ -811,8 +815,20 @@ regrid (bool const do_init)
           // Refinement prolongation must fill all active points
           
           ibset needrecv = box.active + box.overlaps;
+
+          //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) {
+          //    static int center;
+          //    if (!center)
+          //        center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2;
+
+          //    ri->_lower.elt[1] = center;
+          //    ri->_upper.elt[1] = center;
+          //}
           
-          i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
+          int ssize = prolongation_stencil_size(rl);
+          ivect ss = ivect(ssize, 0, ssize);
+          i2vect const stencil_size = i2vect (ss, ss);
+          //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
           
           ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0),
                     "Refinement factors must be integer multiples of each other");
@@ -888,6 +904,15 @@ regrid (bool const do_init)
           // Outer boundaries are synchronised for backward
           // compatibility.
           ibset needrecv = box.ghosts;
+
+          //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) {
+          //    static int center;
+          //    if (!center)
+          //        center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2;
+
+          //    ri->_lower.elt[1] = center;
+          //    ri->_upper.elt[1] = center;
+          //}
 #endif
           
           ibset & sync = box.sync;
@@ -898,14 +923,14 @@ regrid (bool const do_init)
 #if 0
             ibset const ovlp = needrecv & obox.owned;
 #else
-            ibset const ovlp = needrecv & obox.interior;
+            ibset ovlp = needrecv & obox.interior;
 #endif
             
             if (cc == c) {
               ASSERT_cc (ovlp.empty(),
                          "A region may not synchronise from itself");
             }
-            
+
             for (ibset::const_iterator
                    ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
             {
@@ -948,6 +973,7 @@ regrid (bool const do_init)
           // Outer boundaries are synchronised for backward
           // compatibility.
           ibset needrecv = box.ghosts;
+
 #endif
           
           // Points which are synchronised need not be boundary
@@ -960,10 +986,23 @@ regrid (bool const do_init)
           // Prolongation must fill what cannot be synchronised, and
           // also all buffer zones
           needrecv += box.buffers;
+
+          //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri)
+          //{
+          //    static int center;
+          //    if (!center)
+          //        center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2;
+
+          //    ri->_lower.elt[1] = center;
+          //    ri->_upper.elt[1] = center;
+          //}
           
           ibset & bndref = box.bndref;
           
-          i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
+          int ssize = prolongation_stencil_size(rl);
+          ivect ss = ivect(ssize, 0, ssize);
+          i2vect const stencil_size = i2vect (ss, ss);
+          //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
           
           ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0),
                     "Refinement factors must be integer multiples of each other");
@@ -1061,6 +1100,14 @@ regrid (bool const do_init)
             // NOTE: b/c of this we need a low-level sync after the restrict
             needrecv = allrestricted & obox.interior;
           }
+          //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) {
+          //    static int center;
+          //    if (!center)
+          //        center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2;
+
+          //    ri->_lower.elt[1] = center;
+          //    ri->_upper.elt[1] = center;
+          //}
           // Cannot restrict into buffer zones
           assert ((allrestricted & obox.buffers).empty());
           
@@ -1787,7 +1834,10 @@ regrid (bool const do_init)
             // of the new grid structure.  It must fill what cannot be
             // synchronised.
             
-            i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
+            int ssize = prolongation_stencil_size(rl);
+            ivect ss = ivect(ssize, 0, ssize);
+            i2vect const stencil_size = i2vect (ss, ss);
+            //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
             
             ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0),
                       "Refinement factors must be integer multiples of each other");
diff --git a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
index 0dc18ccf6..465dcb5b1 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
@@ -550,18 +550,18 @@ namespace CarpetLib {
     
     
     
-    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 (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");
+    //}
     
     
     
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index c7b2c7581..01d0b95d1 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -44,10 +44,10 @@ class vect {
   
   // Fields
   
+public:
   /** Vector elements.  */
   T elt[D==0 ? 1 : D];
   
-public:
   
   // Constructors