global: Change the way in which the grid hierarchy is stored

Change the way in which the grid hierarchy is stored.  The new hierarchy is
        map
        mglevel
        reflevel
        component
        timelevel
i.e., mglevel moved from the bottom to almost the top.  This is
because mglevel used to be a true multigrid level, but is now meant to
be a convergence level.

Do not allocate all storage all the time.  Allow storage to be
switched on an off per refinement level (and for a single mglevel,
which prompted the change above).  Handle storage management with
CCTK_{In,De}creaseGroupStorage instead of
CCTK_{En,Dis}ableGroupStorage.

darcs-hash:20050201225827-891bb-eae3b6bd092ae8d6b5e49be84c6f09f0e882933e.gz

10 files changed, 438 insertions, 400 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index fff52c35c..dc3354204 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -79,15 +79,16 @@ void dh::recompose (const bool do_prolongate) {
   }
 }
 
-void dh::allocate_bboxes() {
-  boxes.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    boxes.at(rl).resize(h.components(rl));
-    for (int c=0; c<h.components(rl); ++c) {
-      boxes.at(rl).at(c).resize(h.mglevels(rl,c));
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        const ibbox intr = h.extents().at(rl).at(c).at(ml);
-        dboxes & b = boxes.at(rl).at(c).at(ml);
+void dh::allocate_bboxes ()
+{
+  boxes.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    boxes.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      boxes.at(ml).at(rl).resize(h.components(rl));
+      for (int c=0; c<h.components(rl); ++c) {
+        const ibbox intr = h.extents().at(ml).at(rl).at(c);
+        dboxes & b = boxes.at(ml).at(rl).at(c);
 
         // Interior
         // (the interior of the grid has the extent as specified by
@@ -109,21 +110,21 @@ void dh::allocate_bboxes() {
         // (interior + boundaries = exterior)
         b.boundaries = b.exterior - intr;
 
-      } // for ml
-    } // for c
-  } // for rl
+      } // for c
+    } // for rl
+  } // for ml
 }
 
-// Loops over each refinement level, each component, and each multigrid
-// level, executing the "boxesop" member function argument on the corresponding
-// element of the "boxes" member 
-void dh::foreach_reflevel_component_mglevel (dh::boxesop op) {
-
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        dboxes & box = boxes.at(rl).at(c).at(ml);
-        (this->*op)( box, rl, c, ml ); // evaluate member function
+// Loops over each multigrid level, each refinement level, and each
+// component, executing the "boxesop" member function argument on the
+// corresponding element of the "boxes" member
+void dh::foreach_reflevel_component_mglevel (dh::boxesop op)
+{
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      for (int c=0; c<h.components(rl); ++c) {
+        dboxes & box = boxes.at(ml).at(rl).at(c);
+        (this->*op)(box, rl, c, ml); // evaluate member function
       }
     }
   }
@@ -157,8 +158,7 @@ void dh::intersect_sync_with_interior( dh::dboxes & box, int rl, int c, int ml )
 
   // Sync boxes
   for (int cc=0; cc<h.components(rl); ++cc) {
-    dboxes & box1 = boxes.at(rl).at(cc).at(ml);
-    assert (ml<h.mglevels(rl,cc));
+    dboxes & box1 = boxes.at(ml).at(rl).at(cc);
     // intersect boundaries with interior of that component
     ibset ovlp = bnds & box1.interior;
     ovlp.normalize();
@@ -176,7 +176,7 @@ void dh::setup_multigrid_boxes( dh::dboxes & box, int rl, int c, int ml )
 
   // Multigrid boxes
   if (ml>0) {
-    dboxes & bbox = boxes.at(rl).at(c).at(ml-1);
+    dboxes & bbox = boxes.at(ml-1).at(rl).at(c);
     const ibbox intrf = bbox.interior;
     const ibbox extrf = bbox.exterior;
     // Restriction (interior)
@@ -217,7 +217,7 @@ void dh::setup_refinement_interior_boxes( dh::dboxes & box, int rl, int c, int m
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       // Prolongation (interior)
       // TODO: prefer boxes from the same processor
@@ -257,7 +257,7 @@ void dh::setup_refinement_exterior_boxes( dh::dboxes & box, int rl, int c, int m
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       const ibbox& extrf = box1.exterior;
       const ibset& bndsf = box1.boundaries;
@@ -332,7 +332,7 @@ void dh::setup_restrict_interior_boxes( dh::dboxes & box, int rl, int c, int ml
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       // Restriction (interior)
       {
@@ -422,7 +422,7 @@ void dh::assert_assert_assert( dh::dboxes & box, int rl, int c, int ml )
     if (rl==0) {
       const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
       assert (recv_ref_coarse.empty());
-    } else {              // rl!=0
+    } else {                    // rl!=0
       const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
       ibset intr = box.interior;
       ibset received;
@@ -504,23 +504,19 @@ void dh::assert_assert_assert( dh::dboxes & box, int rl, int c, int ml )
 
 void dh::calculate_bases () {
   // Calculate bases
-  bases.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    if (h.components(rl)==0) {
-      bases.at(rl).resize(0);
-    } else {
-      bases.at(rl).resize(h.mglevels(rl,0));
-      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-        dbases & base = bases.at(rl).at(ml);
-        base.exterior = ibbox();
-        base.interior = ibbox();
-        for (int c=0; c<h.components(rl); ++c) {
-          dboxes & box = boxes.at(rl).at(c).at(ml);
-          base.exterior = (base.exterior.expanded_containing(box.exterior));
-          base.interior = (base.interior.expanded_containing(box.interior));
-        }
-        base.boundaries = base.exterior - base.interior;
+  bases.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    bases.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      dbases & base = bases.at(ml).at(rl);
+      base.exterior = ibbox();
+      base.interior = ibbox();
+      for (int c=0; c<h.components(rl); ++c) {
+        dboxes & box = boxes.at(ml).at(rl).at(c);
+        base.exterior = (base.exterior.expanded_containing(box.exterior));
+        base.interior = (base.interior.expanded_containing(box.interior));
       }
+      base.boundaries = base.exterior - base.interior;
     }
   }
 }
@@ -569,8 +565,8 @@ void dh::save_memory ( bool do_prolongate ) {
         assert (! vectorleader);
         (*f)->recompose_free (rl);
       } else {
-        // treat vector groups specially: delete the leader only
-        // after all other elements have been deleted
+        // treat vector groups specially: delete the leader only after
+        // all other elements have been deleted
         if ((*f)->vectorindex == 0) {
           // first element: delete nothing
           if (rl == 0) {
@@ -632,7 +628,7 @@ void dh::do_output_bboxes( dh::dboxes & box, int rl, int c, int ml )
 {
   cout << endl;
   cout << "dh bboxes:" << endl;
-  cout << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
+  cout << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
   cout << "exterior=" << box.exterior << endl;
   cout << "interior=" << box.interior << endl;
   cout << "send_mg_fine=" << box.send_mg_fine << endl;
@@ -652,14 +648,15 @@ void dh::do_output_bboxes( dh::dboxes & box, int rl, int c, int ml )
   cout << "recv_not=" << box.recv_not << endl;
 }
 
-void dh::output_bases () {
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    if (h.components(rl)>0) {
-      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-        dbases & base = bases.at(rl).at(ml);
+void dh::output_bases ()
+{
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      if (h.components(rl)>0) {
+        dbases & base = bases.at(ml).at(rl);
         cout << endl;
         cout << "dh bases:" << endl;
-        cout << "rl=" << rl << " ml=" << ml << endl;
+        cout << "ml=" << ml << " rl=" << rl << endl;
         cout << "exterior=" << base.exterior << endl;
         cout << "interior=" << base.interior << endl;
         cout << "boundaries=" << base.boundaries << endl;
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index a0e8b8e0d..8fa11f152 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -55,7 +55,7 @@ public:
     iblistvect send_ref_coarse;
     iblistvect recv_ref_fine;
     iblistvect recv_ref_coarse;
-    iblistvect send_sync;        // send while syncing
+    iblistvect send_sync;       // send while syncing
     iblistvect send_ref_bnd_fine; // sent to finer grids
     
     ibset boundaries;           // boundaries
@@ -63,14 +63,6 @@ public:
     iblistvect recv_ref_bnd_coarse; // received from coarser grids
     ibset sync_not;             // not received while syncing (outer boundary of that level)
     ibset recv_not;             // not received while syncing or prolongating (globally outer boundary)
-    
-#if 0
-    // after regridding:
-    iblistvect prev_send;       // sent from previous dh
-    iblistvect recv_prev;       // received from previous dh
-    iblistvect send_prev_fine;  // sent to finer
-    iblistvect recv_prev_coarse; // received from coarser
-#endif
   };
   
 private:
@@ -81,19 +73,19 @@ private:
     ibset boundaries;           // boundaries
   };
   
-  typedef vector<dboxes> mboxes; // ... for each multigrid level
-  typedef vector<mboxes> cboxes; // ... for each component
+  typedef vector<dboxes> cboxes; // ... for each component
   typedef vector<cboxes> rboxes; // ... for each refinement level
+  typedef vector<rboxes> mboxes; // ... for each multigrid level
   
-  typedef vector<dbases> mbases; // ... for each multigrid level
-  typedef vector<mbases> rbases; // ... for each refinement level
+  typedef vector<dbases> rbases; // ... for each refinement level
+  typedef vector<rbases> mbases; // ... for each multigrid level
  
-  void allocate_bboxes();
-                                 // generic member function taking a dboxes,
-                                 // a refinement level, a component, and a 
-                                 // multigrid level
-  typedef void    (dh::*boxesop)( dboxes &, int rl, int c, int ml ); 
-  void foreach_reflevel_component_mglevel ( boxesop op );
+  void allocate_bboxes ();
+
+  // generic member function taking a dboxes, a refinement level, a
+  // component, and a multigrid level
+  typedef void (dh::*boxesop) (dboxes &, int rl, int c, int ml); 
+  void foreach_reflevel_component_mglevel (boxesop op);
 
                                  // these all of form 'boxesop'
   void setup_sync_and_refine_boxes( dboxes & b, int rl, int c, int ml );
@@ -120,8 +112,8 @@ public:                         // should be readonly
   int prolongation_order_space; // order of spatial prolongation operator
   int buffer_width;             // buffer inside refined grids
   
-  rboxes boxes;
-  rbases bases;
+  mboxes boxes;
+  mbases bases;
   
   list<ggf*> gfs;               // list of all grid functions
   
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index e493f9491..e3fde8dd0 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -14,12 +14,11 @@ using namespace std;
 template<typename T>
 gf<T>::gf (const int varindex_, const operator_type transport_operator_,
            th& t_, dh& d_,
-           const int tmin_, const int tmax_,
            const int prolongation_order_time_,
            const int vectorlength_, const int vectorindex_,
            gf* const vectorleader_)
   : ggf(varindex_, transport_operator_,
-        t_, d_, tmin_, tmax_, prolongation_order_time_,
+        t_, d_, prolongation_order_time_,
         vectorlength_, vectorindex_, vectorleader_)
 {
   // recompose ();
@@ -51,21 +50,23 @@ gf<T>::~gf () { }
 
 // Access to the data
 template<typename T>
-const data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) const {
-  assert (tl>=tmin && tl<=tmax);
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  return (const data<T>*)storage.at(tl-tmin).at(rl).at(c).at(ml);
+const data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) const
+{
+  assert (tl>=0 and tl<timelevels());
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  return (const data<T>*)storage.at(ml).at(rl).at(c).at(tl);
 }
 
 template<typename T>
-data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) {
-  assert (tl>=tmin && tl<=tmax);
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  return (data<T>*)storage.at(tl-tmin).at(rl).at(c).at(ml);
+data<T>* gf<T>::operator() (int tl, int rl, int c, int ml)
+{
+  assert (tl>=0 and tl<timelevels());
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  return (data<T>*)storage.at(ml).at(rl).at(c).at(tl);
 }
 
 
@@ -76,7 +77,7 @@ ostream& gf<T>::output (ostream& os) const {
   T Tdummy;
   os << "gf<" << typestring(Tdummy) << ">:"
      << varindex << "[" << CCTK_VarName(varindex) << "],"
-     << "dt=[" << tmin << ":" << tmax<< "]";
+     << "tls=" << timelevels();
   return os;
 }
 
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index f3436d43b..fcfb82b3e 100644
--- a/Carpet/CarpetLib/src/gf.hh
+++ b/Carpet/CarpetLib/src/gf.hh
@@ -35,7 +35,7 @@ public:
   // Constructors
   gf (const int varindex, const operator_type transport_operator,
       th& t, dh& d,
-      const int tmin, const int tmax, const int prolongation_order_time,
+      const int prolongation_order_time,
       const int vectorlength, const int vectorindex,
       gf* const vectorleader);
   
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index c778353b8..7387dd298 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -19,24 +19,23 @@ using namespace std;
 // Constructors
 ggf::ggf (const int varindex_, const operator_type transport_operator_,
           th& t_, dh& d_,
-          const int tmin_, const int tmax_,
           const int prolongation_order_time_,
           const int vectorlength_, const int vectorindex_,
           ggf* const vectorleader_)
   : varindex(varindex_), transport_operator(transport_operator_), t(t_),
-    tmin(tmin_), tmax(tmax_),
     prolongation_order_time(prolongation_order_time_),
     h(d_.h), d(d_),
-    storage(tmax-tmin+1),
+    timelevels_(0),
+    storage(h.mglevels()),
     vectorlength(vectorlength_), vectorindex(vectorindex_),
     vectorleader(vectorleader_)
 {
   assert (&t.h == &d.h);
   
   assert (vectorlength >= 1);
-  assert (vectorindex >= 0 && vectorindex < vectorlength);
-  assert ((vectorindex==0 && !vectorleader)
-          || (vectorindex!=0 && vectorleader));
+  assert (vectorindex >= 0 and vectorindex < vectorlength);
+  assert ((vectorindex==0 and !vectorleader)
+          or (vectorindex!=0 and vectorleader));
   
   d.add(this);
 }
@@ -54,20 +53,53 @@ bool ggf::operator== (const ggf& f) const {
 
 
 // Modifiers
+void ggf::set_timelevels (const int ml, const int rl, const int new_timelevels)
+{
+  assert (ml>=0 and ml<(int)storage.size());
+  assert (rl>=0 and rl<(int)storage.at(ml).size());
+  
+  assert (new_timelevels >= 0);
+  
+  if (new_timelevels < timelevels()) {
+    
+    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+      for (int tl=new_timelevels; tl<timelevels(); ++tl) {
+        delete storage.at(ml).at(rl).at(c).at(tl);
+      }
+      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+    } // for c
+    
+  } else if (new_timelevels > timelevels()) {
+    
+    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+      for (int tl=timelevels(); tl<new_timelevels; ++tl) {
+        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
+        storage.at(ml).at(rl).at(c).at(tl)->allocate
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+      } // for tl
+    } // for c
+    
+  }
+  
+  timelevels_ = new_timelevels;
+}
+
+
+
 void ggf::recompose_crop ()
 {
   // Free storage that will not be needed
-  storage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int rl=h.reflevels(); rl<(int)storage.at(tl-tmin).size(); ++rl) {
-      for (int c=0; c<(int)storage.at(tl-tmin).at(rl).size(); ++c) {
-        for (int ml=0; ml<(int)storage.at(tl-tmin).at(rl).at(c).size(); ++ml) {
-          delete storage.at(tl-tmin).at(rl).at(c).at(ml);
-        } // for ml
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=h.reflevels(); rl<(int)storage.at(ml).size(); ++rl) {
+      for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+        for (int tl=0; tl<(int)storage.at(ml).at(rl).at(c).size(); ++tl) {
+          delete storage.at(ml).at(rl).at(c).at(tl);
+        } // for tl
       } // for c
     } // for rl
-    storage.at(tl-tmin).resize(h.reflevels());
-  } // for tl
+    storage.at(ml).resize(h.reflevels());
+  } // for mrl
 }
 
 void ggf::recompose_allocate (const int rl)
@@ -75,58 +107,55 @@ void ggf::recompose_allocate (const int rl)
   // TODO: restructure storage only when needed
   
   // Retain storage that might be needed
-  oldstorage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    oldstorage.at(tl-tmin).resize(h.reflevels());
-    assert (oldstorage.at(tl-tmin).at(rl).size() == 0);
-    oldstorage.at(tl-tmin).at(rl) = storage.at(tl-tmin).at(rl);
-    storage.at(tl-tmin).at(rl).resize(0);
+  oldstorage.resize(storage.size());
+  for (int ml=0; ml<(int)storage.size(); ++ml) {
+    oldstorage.at(ml).resize(storage.at(ml).size());
+    oldstorage.at(ml).at(rl) = storage.at(ml).at(rl);
+    storage.at(ml).at(rl).resize(0);
   }
   
   // Resize structure and allocate storage
-  storage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    storage.at(tl-tmin).resize(h.reflevels());
-    storage.at(tl-tmin).at(rl).resize(h.components(rl));
+  storage.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    storage.at(ml).resize(h.reflevels());
+    storage.at(ml).at(rl).resize(h.components(rl));
     for (int c=0; c<h.components(rl); ++c) {
-      storage.at(tl-tmin).at(rl).at(c).resize(h.mglevels(rl,c));
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        storage.at(tl-tmin).at(rl).at(c).at(ml) = typed_data(tl,rl,c,ml);
-        storage.at(tl-tmin).at(rl).at(c).at(ml)->allocate
-          (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c));
-      } // for ml
+      storage.at(ml).at(rl).at(c).resize(timelevels());
+      for (int tl=0; tl<timelevels(); ++tl) {
+        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
+        storage.at(ml).at(rl).at(c).at(tl)->allocate
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+      } // for tl
     } // for c
-  } // for tl
+  } // for ml
 }
 
 void ggf::recompose_fill (comm_state& state, const int rl,
                           const bool do_prolongate)
 {
   // Initialise the new storage
-  for (int c=0; c<h.components(rl); ++c) {
-    for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      for (int tl=tmin; tl<=tmax; ++tl) {
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int c=0; c<h.components(rl); ++c) {
+      for (int tl=0; tl<timelevels(); ++tl) {
         
         // Find out which regions need to be prolongated
         // (Copy the exterior because some variables are not prolongated)
         // TODO: do this once in the dh instead of for each variable here
-        ibset work (d.boxes.at(rl).at(c).at(ml).exterior);
+        ibset work (d.boxes.at(ml).at(rl).at(c).exterior);
         
         // Copy from old storage, if possible
         // TODO: copy only from interior regions?
-        if (rl<(int)oldstorage.at(tl-tmin).size()) {
-          for (int cc=0; cc<(int)oldstorage.at(tl-tmin).at(rl).size(); ++cc) {
-            if (ml<(int)oldstorage.at(tl-tmin).at(rl).at(cc).size()) {
-              // TODO: prefer same processor, etc., see dh.cc
-              ibset ovlp
-                = work & oldstorage.at(tl-tmin).at(rl).at(cc).at(ml)->extent();
-              ovlp.normalize();
-              work -= ovlp;
-              for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
-                storage.at(tl-tmin).at(rl).at(c).at(ml)->copy_from
-                  (state, oldstorage.at(tl-tmin).at(rl).at(cc).at(ml), *r);
-              }
-            } // if ml
+        if (rl<(int)oldstorage.at(ml).size()) {
+          for (int cc=0; cc<(int)oldstorage.at(ml).at(rl).size(); ++cc) {
+            // TODO: prefer same processor, etc., see dh.cc
+            ibset ovlp
+              = work & oldstorage.at(ml).at(rl).at(cc).at(tl)->extent();
+            ovlp.normalize();
+            work -= ovlp;
+            for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
+              storage.at(ml).at(rl).at(c).at(tl)->copy_from
+                (state, oldstorage.at(ml).at(rl).at(cc).at(tl), *r);
+            }
           } // for cc
         } // if rl
         
@@ -135,18 +164,17 @@ void ggf::recompose_fill (comm_state& state, const int rl,
           if (rl>0) {
             if (transport_operator != op_none) {
               const int numtl = prolongation_order_time+1;
-              assert (tmax-tmin+1 >= numtl);
+              assert (timelevels() >= numtl);
               vector<int> tls(numtl);
               vector<CCTK_REAL> times(numtl);
               for (int i=0; i<numtl; ++i) {
-                tls.at(i) = tmax - i;
+                tls.at(i) = i;
                 times.at(i) = t.time(tls.at(i),rl-1,ml);
               }
-              for (int cc=0; cc<(int)storage.at(tl-tmin).at(rl-1).size(); ++cc) {
+              for (int cc=0; cc<(int)storage.at(ml).at(rl-1).size(); ++cc) {
                 vector<const gdata*> gsrcs(numtl);
                 for (int i=0; i<numtl; ++i) {
-                  gsrcs.at(i)
-                    = storage.at(tls.at(i)-tmin).at(rl-1).at(cc).at(ml);
+                  gsrcs.at(i) = storage.at(ml).at(rl-1).at(cc).at(tls.at(i));
                   assert (gsrcs.at(i)->extent() == gsrcs.at(0)->extent());
                 }
                 const CCTK_REAL time = t.time(tl,rl,ml);
@@ -154,13 +182,17 @@ void ggf::recompose_fill (comm_state& state, const int rl,
                 // TODO: choose larger regions first
                 // TODO: prefer regions from the same processor
                 const iblist& list
-                  = d.boxes.at(rl).at(c).at(ml).recv_ref_coarse.at(cc);
-                for (iblist::const_iterator iter=list.begin(); iter!=list.end(); ++iter) {
+                  = d.boxes.at(ml).at(rl).at(c).recv_ref_coarse.at(cc);
+                for (iblist::const_iterator iter=list.begin();
+                     iter!=list.end(); ++iter)
+                {
                   ibset ovlp = work & *iter;
                   ovlp.normalize();
                   work -= ovlp;
-                  for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
-                    storage.at(tl-tmin).at(rl).at(c).at(ml)->interpolate_from
+                  for (ibset::const_iterator r=ovlp.begin();
+                       r!=ovlp.end(); ++r)
+                  {
+                    storage.at(ml).at(rl).at(c).at(tl)->interpolate_from
                       (state, gsrcs, times, *r, time,
                        d.prolongation_order_space, prolongation_order_time);
                   } // for r
@@ -176,58 +208,58 @@ void ggf::recompose_fill (comm_state& state, const int rl,
         // TODO: check this.
         
       } // for tl
-    } // for ml
-  } // for c
+    } // for c
+  } // for ml
 }
 
 void ggf::recompose_free (const int rl)
 {
   // Delete old storage
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int c=0; c<(int)oldstorage.at(tl-tmin).at(rl).size(); ++c) {
-      for (int ml=0; ml<(int)oldstorage.at(tl-tmin).at(rl).at(c).size(); ++ml) {
-        delete oldstorage.at(tl-tmin).at(rl).at(c).at(ml);
-      } // for ml
+  for (int ml=0; ml<(int)oldstorage.size(); ++ml) {
+    for (int c=0; c<(int)oldstorage.at(ml).at(rl).size(); ++c) {
+      for (int tl=0; tl<timelevels(); ++tl) {
+        delete oldstorage.at(ml).at(rl).at(c).at(tl);
+      } // for tl
     } // for c
-    oldstorage.at(tl-tmin).at(rl).resize(0);
-  } // for tl
+    oldstorage.at(ml).at(rl).resize(0);
+  } // for ml
 }
 
 void ggf::recompose_bnd_prolongate (comm_state& state, const int rl,
-                                       const bool do_prolongate)
+                                    const bool do_prolongate)
 {
   if (do_prolongate) {
     // Set boundaries
     if (rl>0) {
-      for (int c=0; c<h.components(rl); ++c) {
-        for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-          for (int tl=tmin; tl<=tmax; ++tl) {
+      for (int ml=0; ml<h.mglevels(); ++ml) {
+        for (int c=0; c<h.components(rl); ++c) {
+          for (int tl=0; tl<timelevels(); ++tl) {
             
             // TODO: assert that reflevel 0 boundaries are copied
             const CCTK_REAL time = t.time(tl,rl,ml);
             ref_bnd_prolongate (state,tl,rl,c,ml,time);
             
           } // for tl
-        } // for ml
-      } // for c
+        } // for c
+      } // for ml
     } // if rl
   } // if do_prolongate
 }
 
 void ggf::recompose_sync (comm_state& state, const int rl,
-                             const bool do_prolongate)
+                          const bool do_prolongate)
 {
   if (do_prolongate) {
     // Set boundaries
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        for (int tl=tmin; tl<=tmax; ++tl) {
+    for (int ml=0; ml<h.mglevels(); ++ml) {
+      for (int c=0; c<h.components(rl); ++c) {
+        for (int tl=0; tl<timelevels(); ++tl) {
           
           sync (state,tl,rl,c,ml);
           
         } // for tl
-      } // for ml
-    } // for c
+      } // for c
+    } // for ml
   } // if do_prolongate
 }
 
@@ -235,28 +267,28 @@ void ggf::recompose_sync (comm_state& state, const int rl,
 
 // Cycle the time levels by rotating the data sets
 void ggf::cycle (int rl, int c, int ml) {
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  gdata* tmpdata = storage.at(tmin-tmin).at(rl).at(c).at(ml);
-  for (int tl=tmin; tl<=tmax-1; ++tl) {
-    storage.at(tl-tmin).at(rl).at(c).at(ml) = storage.at(tl+1-tmin).at(rl).at(c).at(ml);
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  gdata* tmpdata = storage.at(ml).at(rl).at(c).at(timelevels()-1);
+  for (int tl=timelevels()-1; tl>0; --tl) {
+    storage.at(ml).at(rl).at(c).at(tl) = storage.at(ml).at(rl).at(c).at(tl-1);
   }
-  storage.at(tmax-tmin).at(rl).at(c).at(ml) = tmpdata;
+  storage.at(ml).at(rl).at(c).at(0) = tmpdata;
 }
 
 // Flip the time levels by exchanging the data sets
 void ggf::flip (int rl, int c, int ml) {
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  for (int tl=0; tl<(tmax-tmin)/2; ++tl) {
-    const int tl1 = tmin + tl;
-    const int tl2 = tmax - tl;
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  for (int tl=0; tl<(timelevels()-1)/2; ++tl) {
+    const int tl1 =                  tl;
+    const int tl2 = timelevels()-1 - tl;
     assert (tl1 < tl2);
-    gdata* tmpdata = storage.at(tl1-tmin).at(rl).at(c).at(ml);
-    storage.at(tl1-tmin).at(rl).at(c).at(ml) = storage.at(tl2-tmin).at(rl).at(c).at(ml);
-    storage.at(tl2-tmin).at(rl).at(c).at(ml) = tmpdata;
+    gdata* tmpdata = storage.at(ml).at(rl).at(c).at(tl1);
+    storage.at(ml).at(rl).at(c).at(tl1) = storage.at(ml).at(rl).at(c).at(tl2);
+    storage.at(ml).at(rl).at(c).at(tl2) = tmpdata;
   }
 }
 
@@ -271,21 +303,21 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const ibbox dh::dboxes::* send_box)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
-  const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_box;
-  const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_box;
+  assert (ml2<h.mglevels());
+  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_box;
+  const ibbox send = d.boxes.at(ml2).at(rl2).at(c2).*send_box;
   assert (all(recv.shape()==send.shape()));
   // copy the content
   assert (recv==send);
-  storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-    (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), recv);
+  storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+    (state, storage.at(ml2).at(rl2).at(c2).at(tl2), recv);
 }
 
 // Copy regions
@@ -295,24 +327,24 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const iblist dh::dboxes::* send_list)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
-  const iblist recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const iblist send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  assert (ml2<h.mglevels());
+  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const iblist send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (recv.size()==send.size());
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(), s=send.begin();
        r!=recv.end(); ++r, ++s) {
     // (use the send boxes for communication)
     // copy the content
-    storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-      (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r);
+    storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+      (state, storage.at(ml2).at(rl2).at(c2).at(tl2), *r);
   }
 }
 
@@ -323,25 +355,25 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const iblistvect dh::dboxes::* send_listvect)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
-    assert (ml2<h.mglevels(rl2,c2));
-    const iblist recv = (d.boxes.at(rl1).at(c1).at(ml1).*recv_listvect).at(c2);
-    const iblist send = (d.boxes.at(rl2).at(c2).at(ml2).*send_listvect).at(c1);
+    assert (ml2<h.mglevels());
+    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+    const iblist send = (d.boxes.at(ml2).at(rl2).at(c2).*send_listvect).at(c1);
     assert (recv.size()==send.size());
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(), s=send.begin();
       	 r!=recv.end(); ++r, ++s) {
       // (use the send boxes for communication)
       // copy the content
-      storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-      	(state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r);
+      storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+      	(state, storage.at(ml2).at(rl2).at(c2).at(tl2), *r);
     }
   }
 }
@@ -354,33 +386,30 @@ void ggf::intercat (comm_state& state,
                     const ibbox dh::dboxes::* send_list,
                     CCTK_REAL time)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 and ml2<h.mglevels());
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-    assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-    assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-    gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
     times.at(i) = t.time(tl2s.at(i),rl2,ml2);
   }
   
-  const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const ibbox send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (all(recv.shape()==send.shape()));
   // interpolate the content
   assert (recv==send);
-  storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+  storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
     (state, gsrcs, times, recv, time,
      d.prolongation_order_space, prolongation_order_time);
 }
@@ -393,36 +422,33 @@ void ggf::intercat (comm_state& state,
                     const iblist dh::dboxes::* send_list,
                     const CCTK_REAL time)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 and ml2<h.mglevels());
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-    assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-    assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-    gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
     times.at(i) = t.time(tl2s.at(i),rl2,ml2);
   }
   
-  const iblist recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const iblist send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const iblist send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (recv.size()==send.size());
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(), s=send.begin();
        r!=recv.end(); ++r, ++s) {
     // (use the send boxes for communication)
     // interpolate the content
-    storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+    storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
       (state, gsrcs, times, *r, time,
        d.prolongation_order_space, prolongation_order_time);
   }
@@ -436,37 +462,34 @@ void ggf::intercat (comm_state& state,
                     const iblistvect dh::dboxes::* send_listvect,
                     const CCTK_REAL time)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
-    assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+    assert (ml2>=0 and ml2<h.mglevels());
     
     vector<const gdata*> gsrcs(tl2s.size());
     vector<CCTK_REAL> times(tl2s.size());
     for (int i=0; i<(int)gsrcs.size(); ++i) {
-      assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-      assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-      assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-      gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+      gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
       times.at(i) = t.time(tl2s.at(i),rl2,ml2);
     }
     
-    const iblist recv = (d.boxes.at(rl1).at(c1).at(ml1).*recv_listvect).at(c2);
-    const iblist send = (d.boxes.at(rl2).at(c2).at(ml2).*send_listvect).at(c1);
+    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+    const iblist send = (d.boxes.at(ml2).at(rl2).at(c2).*send_listvect).at(c1);
     assert (recv.size()==send.size());
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(), s=send.begin();
       	 r!=recv.end(); ++r, ++s) {
       // (use the send boxes for communication)
       // interpolate the content
-      storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+      storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
       	(state, gsrcs, times, *r, time,
 	 d.prolongation_order_space, prolongation_order_time);
     }
@@ -503,9 +526,9 @@ void ggf::ref_bnd_prolongate (comm_state& state,
   if (transport_operator == op_none) return;
   vector<int> tl2s;
   // Interpolation in time
-  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  assert (timelevels() >= prolongation_order_time+1);
   tl2s.resize(prolongation_order_time+1);
-  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i;
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse,
 	    tl2s,rl-1,  ml, &dh::dboxes::send_ref_bnd_fine,
@@ -567,9 +590,9 @@ void ggf::ref_prolongate (comm_state& state,
   if (transport_operator == op_none) return;
   vector<int> tl2s;
   // Interpolation in time
-  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  assert (timelevels() >= prolongation_order_time+1);
   tl2s.resize(prolongation_order_time+1);
-  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i;
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_coarse,
 	    tl2s,rl-1,  ml, &dh::dboxes::send_ref_fine,
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index f1ac54baa..0fcaf960f 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -34,10 +34,10 @@ class ggf {
   typedef vector<iblist> iblistvect;
   
   typedef gdata*        tdata;  // data ...
-  typedef vector<tdata> mdata;  // ... for each multigrid level
-  typedef vector<mdata> cdata;  // ... for each component
+  typedef vector<tdata> fdata;  // ... for each time level
+  typedef vector<fdata> cdata;  // ... for each component
   typedef vector<cdata> rdata;  // ... for each refinement level
-  typedef vector<rdata> fdata;  // ... for each time level
+  typedef vector<rdata> mdata;  // ... for each multigrid level
   
 public:				// should be readonly
   
@@ -45,15 +45,17 @@ public:				// should be readonly
   const int varindex;           // Cactus variable index
   const operator_type transport_operator;
   
-  const th &t;               // time hierarchy
-  const int tmin, tmax;		// timelevels
+  const th &t;                  // time hierarchy
   const int prolongation_order_time; // order of temporal prolongation operator
   
-  const gh &h;               // grid hierarchy
+  const gh &h;                  // grid hierarchy
   dh &d;			// data hierarchy
 
+private:
+  int timelevels_;              // time levels
+
 protected:
-  fdata storage;		// storage
+  mdata storage;		// storage
   
 public:
   const int vectorlength;       // vector length
@@ -61,14 +63,13 @@ public:
   const ggf* vectorleader;      // first vector element
   
 private:
-  fdata oldstorage;
+  mdata oldstorage;            // temporary storage
   
 public:
 
   // Constructors
   ggf (const int varindex, const operator_type transport_operator,
        th& t, dh& d,
-       const int tmin, const int tmax,
        const int prolongation_order_time,
        const int vectorlength, const int vectorindex,
        ggf* const vectorleader);
@@ -78,11 +79,18 @@ public:
 
   // Comparison
   bool operator== (const ggf& f) const;
-
-
-
+  
+  // Querying
+  int timelevels () const
+  {
+    return timelevels_;
+  }
+  
+  
+  
   // Modifiers
-  // void recompose ();
+  void set_timelevels (int ml, int rl, int new_timelevels);
+
   void recompose_crop ();
   void recompose_allocate (int rl);
   void recompose_fill (comm_state& state, int rl, bool do_prolongate);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index c73bc9dd3..7b5f97c81 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -29,7 +29,7 @@ gh::gh (const int reffact_, const centering refcent_,
 gh::~gh () { }
 
 // Modifiers
-void gh::recompose (const rexts& exts,
+void gh::recompose (const mexts& exts,
                     const rbnds& outer_bounds,
                     const rprocs& procs,
                     const bool do_prolongate)
@@ -70,45 +70,47 @@ void gh::recompose (const rexts& exts,
 
 void gh::check_processor_number_consistency () {
   for (int rl=0; rl<reflevels(); ++rl) {
-    assert (processors().size() == extents().size());
-    assert (outer_boundaries().size() == extents().size());
+    assert (processors().size() == extents().at(0).size());
+    assert (outer_boundaries().size() == extents().at(0).size());
     for (int c=0; c<components(rl); ++c) {
-      assert (processors().at(rl).size() == extents().at(rl).size());
-      assert (outer_boundaries().at(rl).size() == extents().at(rl).size());
+      assert (processors().at(rl).size() == extents().at(0).at(rl).size());
+      assert (outer_boundaries().at(rl).size() == extents().at(0).at(rl).size());
     }
   }
 }
   
-void gh::check_multigrid_consistency () {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    for (int c=0; c<components(rl); ++c) {
-      assert (mglevels(rl,c)>0);
-      for (int ml=1; ml<mglevels(rl,c); ++ml) {
-	assert (all(extents().at(rl).at(c).at(ml).stride()
-		    == ivect(mgfact) * extents().at(rl).at(c).at(ml-1).stride()));
+void gh::check_multigrid_consistency ()
+{
+  assert (mglevels()>0);
+  for (int ml=1; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      for (int c=0; c<components(rl); ++c) {
+	assert (all(extents().at(ml).at(rl).at(c).stride()
+		    == ivect(mgfact) * extents().at(ml-1).at(rl).at(c).stride()));
         // TODO: put the check back in, taking outer boundaries into
         // account
 #if 0
-	assert (extents().at(rl).at(c).at(ml)
-		.contracted_for(extents().at(rl).at(c).at(ml-1))
-		.is_contained_in(extents().at(rl).at(c).at(ml-1)));
+	assert (extents().at(ml).at(rl).at(c)
+		.contracted_for(extents().at(ml-1).at(rl).at(c))
+		.is_contained_in(extents().at(ml-1).at(rl).at(c)));
 #endif
       }
     }
   }
 }
   
-void gh::check_component_consistency () {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    assert (components(rl)>0);
-    for (int c=0; c<components(rl); ++c) {
-      for (int ml=0; ml<mglevels(rl,c); ++ml) {
-        const ibbox &b  = extents().at(rl).at(c).at(ml);
-        const ibbox &b0 = extents().at(rl).at(0).at(ml);
+void gh::check_component_consistency ()
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      assert (components(rl)>0);
+      for (int c=0; c<components(rl); ++c) {
+        const ibbox &b  = extents().at(ml).at(rl).at(c);
+        const ibbox &b0 = extents().at(ml).at(rl).at(0);
 	assert (all(b.stride() == b0.stride()));
 	assert (b.is_aligned_with(b0));
         for (int cc=c+1; cc<components(rl); ++cc) {
-          assert ((b & extents().at(rl).at(cc).at(ml)).empty());
+          assert ((b & extents().at(ml).at(rl).at(cc)).empty());
         }
       }
     }
@@ -127,40 +129,40 @@ void gh::check_base_grid_extent () {
   }
 }
   
-void gh::check_refinement_levels () {
-  for (int rl=1; rl<reflevels(); ++rl) {
-    assert (all(extents().at(rl-1).at(0).at(0).stride()
-		== ivect(reffact) * extents().at(rl).at(0).at(0).stride()));
-    // Check contained-ness:
-    // first take all coarse grids ...
-    ibset all;
-    for (int c=0; c<components(rl-1); ++c) {
-      all |= extents().at(rl-1).at(c).at(0);
-    }
-    // ... remember their size ...
-    const int sz = all.size();
-    // ... then add the coarsified fine grids ...
-    for (int c=0; c<components(rl); ++c) {
-      all |= extents().at(rl).at(c).at(0).contracted_for(extents().at(rl-1).at(0).at(0));
+void gh::check_refinement_levels ()
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=1; rl<reflevels(); ++rl) {
+      assert (all(extents().at(ml).at(rl-1).at(0).stride()
+                  == ivect(reffact) * extents().at(ml).at(rl).at(0).stride()));
+      // Check contained-ness:
+      // first take all coarse grids ...
+      ibset all;
+      for (int c=0; c<components(rl-1); ++c) {
+        all |= extents().at(ml).at(rl-1).at(c);
+      }
+      // ... remember their size ...
+      const int sz = all.size();
+      // ... then add the coarsified fine grids ...
+      for (int c=0; c<components(rl); ++c) {
+        all |= extents().at(ml).at(rl).at(c).contracted_for(extents().at(ml).at(rl-1).at(0));
+      }
+      // ... and then check the sizes:
+      assert (all.size() == sz);
     }
-    // ... and then check the sizes:
-    assert (all.size() == sz);
   }
 }
 
-void gh::calculate_base_extents_of_all_levels () {
-  _bases.resize(reflevels());
-  for (int rl=0; rl<reflevels(); ++rl) {
-    if (components(rl)==0) {
-      _bases.at(rl).resize(0);
-    } else {
-      _bases.at(rl).resize(mglevels(rl,0));
-      for (int ml=0; ml<mglevels(rl,0); ++ml) {
-        _bases.at(rl).at(ml) = ibbox();
-        ibbox &bb = _bases.at(rl).at(ml);
-	for (int c=0; c<components(rl); ++c) {
-	  bb = bb.expanded_containing(extents().at(rl).at(c).at(ml));
-	}
+void gh::calculate_base_extents_of_all_levels ()
+{
+  _bases.resize(mglevels());
+  for (int ml=0; ml<mglevels(); ++ml) {
+    _bases.at(ml).resize(reflevels());
+    for (int rl=0; rl<reflevels(); ++rl) {
+      _bases.at(ml).at(rl) = ibbox();
+      ibbox &bb = _bases.at(ml).at(rl);
+      for (int c=0; c<components(rl); ++c) {
+        bb = bb.expanded_containing(extents().at(ml).at(rl).at(c));
       }
     }
   }
@@ -198,14 +200,15 @@ void gh::remove (dh* d) {
 }
 
 
-void gh::do_output_bboxes (ostream& os) const {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    for (int c=0; c<components(rl); ++c) {
-      for (int ml=0; ml<mglevels(rl,c); ++ml) {
+void gh::do_output_bboxes (ostream& os) const
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      for (int c=0; c<components(rl); ++c) {
         os << endl;
         os << "gh bboxes:" << endl;
-        os << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
-        os << "extent=" << extents().at(rl).at(c).at(ml) << endl;
+        os << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
+        os << "extent=" << extents().at(ml).at(rl).at(c) << endl;
         os << "outer_boundary=" << outer_boundaries().at(rl).at(c) << endl;
         os << "processor=" << processors().at(rl).at(c) << endl;
       }
@@ -213,14 +216,15 @@ void gh::do_output_bboxes (ostream& os) const {
   }
 }
 
-void gh::do_output_bases (ostream& os) const {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    if (components(rl)>0) {
-      for (int ml=0; ml<mglevels(rl,0); ++ml) {
+void gh::do_output_bases (ostream& os) const
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      if (components(rl)>0) {
         os << endl;
         os << "gh bases:" << endl;
-        os << "rl=" << rl << " ml=" << ml << endl;
-        os << "base=" << bases().at(rl).at(ml) << endl;
+        os << "ml=" << ml << " rl=" << rl << endl;
+        os << "base=" << bases().at(ml).at(rl) << endl;
       }
     }
   }
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 6a7bb2e83..3b7f0a2c8 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -33,9 +33,9 @@ class gh {
 public:
   
   // Types
-  typedef vector<ibbox> mexts;	// ... for each multigrid level
-  typedef vector<mexts> cexts;	// ... for each component
+  typedef vector<ibbox> cexts;	// ... for each component
   typedef vector<cexts> rexts;	// ... for each refinement level
+  typedef vector<rexts> mexts;	// ... for each multigrid level
   
   typedef vector<bbvect> cbnds;	// ... for each component
   typedef vector<cbnds> rbnds;	// ... for each refinement level
@@ -56,9 +56,9 @@ public:				// should be readonly
   
   
 private:
-  vector<vector<ibbox> > _bases; // [rl][ml]
+  vector<vector<ibbox> > _bases; // [ml][rl]
   // TODO: invent structure for this
-  rexts _extents;		// extents of all grids
+  mexts _extents;		// extents of all grids
   rbnds _outer_boundaries;	// boundary descriptions of all grids
   rprocs _processors;		// processor numbers of all grids
 
@@ -76,11 +76,12 @@ public:
   virtual ~gh ();
   
   // Modifiers
-  void recompose (const rexts& exts, const rbnds& outer_bounds,
+  void recompose (const mexts& exts, const rbnds& outer_bounds,
 		  const rprocs& procs,
                   const bool do_prolongate);
 
-  const rexts & extents() const {
+  const mexts & extents() const
+  {
     return _extents;
   }
 
@@ -97,16 +98,20 @@ public:
   }
   
   // Accessors
-  int reflevels () const {
+  int mglevels () const
+  {
     return (int)_extents.size();
   }
   
-  int components (const int rl) const {
-    return (int)_extents.at(rl).size();
+  int reflevels () const
+  {
+    if (mglevels() == 0) return 0;
+    return (int)_extents.at(0).size();
   }
   
-  int mglevels (const int rl, const int c) const {
-    return (int)_extents.at(rl).at(c).size();
+  int components (const int rl) const
+  {
+    return (int)_extents.at(0).at(rl).size();
   }
   
   bbvect outer_boundary (const int rl, const int c) const {
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index 18559ae75..77415fffc 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -25,29 +25,29 @@ th::~th () {
 }
 
 // Modifiers
-void th::recompose () {
-  times.resize(h.reflevels());
-  deltas.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    const int old_mglevels = times.at(rl).size();
-    CCTK_REAL mgtime;
-    // Select default time
-    if (old_mglevels==0 && rl==0) {
-      mgtime = 0;
-    } else if (old_mglevels==0) {
-      mgtime = times.at(rl-1).at(0);
-    } else {
-      mgtime = times.at(rl).at(old_mglevels-1);
-    }
-    times.at(rl).resize(h.mglevels(rl,0), mgtime);
-    deltas.at(rl).resize(h.mglevels(rl,0));
-    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-      if (rl==0 && ml==0) {
-	deltas.at(rl).at(ml) = delta;
+void th::recompose ()
+{
+  const int old_mglevels = times.size();
+  times.resize(h.mglevels());
+  deltas.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    const int old_reflevels = times.at(ml).size();
+    times.at(ml).resize(h.reflevels());
+    deltas.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      if (old_mglevels==0) {
+        times.at(ml).at(rl) = 0;
+      } else if (rl < old_reflevels) {
+        // do nothing
+      } else {
+        times.at(ml).at(rl) = times.at(ml).at(rl-1);
+      }
+      if (ml==0 and rl==0) {
+	deltas.at(ml).at(rl) = delta;
       } else if (ml==0) {
-	deltas.at(rl).at(ml) = deltas.at(rl-1).at(ml) / h.reffact;
+	deltas.at(ml).at(rl) = deltas.at(ml).at(rl-1) / h.reffact;
       } else {
-	deltas.at(rl).at(ml) = deltas.at(rl).at(ml-1) * h.mgfact;
+	deltas.at(ml).at(rl) = deltas.at(ml-1).at(rl) * h.mgfact;
       }
     }
   }
@@ -59,11 +59,13 @@ void th::recompose () {
 void th::output (ostream& os) const {
   os << "th:"
      << "times={";
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-      if (!(rl==0 && ml==0)) os << ",";
-      os << rl << ":" << ml << ":"
-	 << times.at(rl).at(ml) << "(" << deltas.at(rl).at(ml) << ")";
+  const char * sep = "";
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      os << sep
+         << ml << ":" << rl << ":"
+	 << times.at(ml).at(rl) << "(" << deltas.at(ml).at(rl) << ")";
+      sep = ",";
     }
   }
   os << "}";
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 6ef1308aa..42bd6504e 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -33,8 +33,8 @@ public:				// should be readonly
 private:
   
   CCTK_REAL delta;		// time step
-  vector<vector<CCTK_REAL> > times; // current times
-  vector<vector<CCTK_REAL> > deltas; // time steps
+  vector<vector<CCTK_REAL> > times; // current times [ml][rl]
+  vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
   
 public:
   
@@ -48,38 +48,44 @@ public:
   void recompose ();
   
   // Time management
-  CCTK_REAL get_time (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return times.at(rl).at(ml);
+  CCTK_REAL get_time (const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return times.at(ml).at(rl);
   }
   
-  void set_time (const int rl, const int ml, const CCTK_REAL t) {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    times.at(rl).at(ml) = t;
+  void set_time (const int rl, const int ml, const CCTK_REAL t)
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    times.at(ml).at(rl) = t;
   }
   
-  void advance_time (const int rl, const int ml) {
+  void advance_time (const int rl, const int ml)
+  {
     set_time(rl,ml, get_time(rl,ml) + get_delta(rl,ml));
   }
   
-  CCTK_REAL get_delta (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return deltas.at(rl).at(ml);
+  CCTK_REAL get_delta (const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return deltas.at(ml).at(rl);
   }
   
-  void set_delta (const int rl, const int ml, const CCTK_REAL dt) {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    deltas.at(rl).at(ml) = dt;
+  void set_delta (const int rl, const int ml, const CCTK_REAL dt)
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    deltas.at(ml).at(rl) = dt;
   }
   
-  CCTK_REAL time (const int tl, const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return get_time(rl, ml) + tl * get_delta(rl, ml);
+  CCTK_REAL time (const int tl, const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return get_time(rl, ml) - tl * get_delta(rl, ml);
   }
   
   // Output