CarpetLib: Split recompose functions into regrid and recompose

Split recompose functions into two stages, regrid and recompose.  The
first stage, regrid, changes the grid structure in the gh and dh
classes.  The second stage, recompose, changes the values of the
actual grid functions, i.e., changes the gf<T> and data<T> objects.
The second stage has to be called individually for every refinement
level.

This is necessary since the boundary conditions need to be applied
after recomposing one refinement level, before the next fine
refinement level can be recomposed.

darcs-hash:20060904230433-dae7b-3ba1982460f57b34da11a6fbb6b4b524dc5b348f.gz

9 files changed, 137 insertions, 132 deletions

diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index 54501823e..1a064d0bb 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -16,14 +16,14 @@ BOOLEAN barriers "Insert barriers at strategic places for debugging purposes (sl
 
 
 
-BOOLEAN poison_new_memory "Try to catch uninitialised data by setting newly allocated memory to values that will catch your attention" STEERABLE=always
+BOOLEAN omit_prolongation_points_when_restricting "Do not restrict to points which are used to prolongate the boundary" STEERABLE=always
 {
 } "no"
 
-CCTK_INT poison_value "Integer value (0..255) used to poison new timelevels (with memset)" STEERABLE=always
+INT proper_nesting_distance "Minimum distance (in grid points) between two level interfaces" STEERABLE=always
 {
-  0:255 :: "Must fit into a byte.  Use 0 for zero, 255 for nan, and e.g. 113 for a large value."
-} 255
+  0:* :: "any non-negative value is fine; the default value is just a guess"
+} 4
 
 
 
@@ -31,34 +31,28 @@ BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=
 {
 } "no"
 
-BOOLEAN print_timestats "Print timing statistics at every iteration" STEERABLE=always
-{
-} "no"
 
-INT print_timestats_every "Print timing statistics periodically" STEERABLE=always
-{
-  0   :: "don't report"
-  1:* :: "report every so many iterations"
-} 0
 
-BOOLEAN save_memory_during_regridding "Save some memory during regridding at the expense of speed" STEERABLE=always
+BOOLEAN poison_new_memory "Try to catch uninitialised data by setting newly allocated memory to values that will catch your attention" STEERABLE=always
 {
-} "yes"
+} "no"
 
-BOOLEAN fast_recomposing "Take shortcuts during recomposing (EXPERIMENTAL)" STEERABLE=always
+CCTK_INT poison_value "Integer value (0..255) used to poison new timelevels (with memset)" STEERABLE=always
 {
-} "no"
+  0:255 :: "Must fit into a byte.  Use 0 for zero, 255 for nan, and e.g. 113 for a large value."
+} 255
 
 
 
-BOOLEAN omit_prolongation_points_when_restricting "Do not restrict to points which are used to prolongate the boundary" STEERABLE=always
+BOOLEAN print_timestats "Print timing statistics at every iteration" STEERABLE=always
 {
 } "no"
 
-INT proper_nesting_distance "Minimum distance (in grid points) between two level interfaces" STEERABLE=always
+INT print_timestats_every "Print timing statistics periodically" STEERABLE=always
 {
-  0:* :: "any non-negative value is fine; the default value is just a guess"
-} 4
+  0   :: "don't report"
+  1:* :: "report every so many iterations"
+} 0
 
 
 
@@ -80,6 +74,8 @@ STRING memstat_file "File name in which memstat output is collected (because std
   "^.+$" :: "file name"
 } "carpetlib-memory-statistics"
 
+
+
 SHARES: IO
 
 USES STRING out_dir
@@ -117,3 +113,11 @@ BOOLEAN use_collective_communication_buffers "Use collective buffers for MPI com
 BOOLEAN minimise_outstanding_communications "Minimise the number of Isend/Irecv operations that are submitted concurrently -- DEPRECATED - DO NOT USE ANYMORE" STEERABLE=always
 {
 } "no"
+
+BOOLEAN save_memory_during_regridding "Save some memory during regridding at the expense of speed -- DEPRECATED - DO NOT USE ANYMORE" STEERABLE=always
+{
+} "yes"
+
+BOOLEAN fast_recomposing "Take shortcuts during recomposing -- DEPRECATED - DO NOT USE ANYMORE" STEERABLE=always
+{
+} "no"
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index f43c20efa..a10fa1450 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -31,7 +31,10 @@ dh::dh (gh& h_,
   assert (all(buffers[0]>=0 and buffers[1]>=0));
   h.add(this);
   CHECKPOINT;
-  recompose (false);
+  regrid ();
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    recompose (rl, false);
+  }
 }
 
 // Destructors
@@ -49,12 +52,12 @@ int dh::prolongation_stencil_size () const
 }
 
 // Modifiers
-void dh::recompose (const bool do_prolongate)
+void dh::regrid ()
 {
   DECLARE_CCTK_PARAMETERS;
   
   CHECKPOINT;
-
+  
   boxes.clear();
 
   allocate_bboxes();
@@ -76,12 +79,29 @@ void dh::recompose (const bool do_prolongate)
   }
   
   foreach_reflevel_component_mglevel (&dh::check_bboxes);
+}
 
-  if (not save_memory_during_regridding) {
-    save_time(do_prolongate);
-  } else {
-    save_memory(do_prolongate);
+void dh::recompose (const int rl, const bool do_prolongate)
+{
+  assert (rl>=0 and rl<h.reflevels());
+  
+  for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+    (*f)->recompose_crop ();
   }
+  
+  for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+    (*f)->recompose_allocate (rl);
+    for (comm_state state; not state.done(); state.step()) {
+      (*f)->recompose_fill (state, rl, do_prolongate);
+    }
+    (*f)->recompose_free (rl);
+    for (comm_state state; not state.done(); state.step()) {
+      (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
+    }
+    for (comm_state state; not state.done(); state.step()) {
+      (*f)->recompose_sync (state, rl, do_prolongate);
+    }
+  } // for all grid functions of same vartype
 }
 
 void dh::allocate_bboxes ()
@@ -717,72 +737,7 @@ void dh::calculate_bases ()
   }
 }
 
-void dh::save_time (bool do_prolongate)
-{
-  DECLARE_CCTK_PARAMETERS;
 
-  for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-    (*f)->recompose_crop ();
-  }
-  bool any_level_changed = false;
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    
-    // TODO: calculate this_level_changed only once
-    bool const this_level_changed
-      = gfs.begin()!=gfs.end() and (*gfs.begin())->recompose_did_change (rl);
-    any_level_changed |= this_level_changed;
-    if (not fast_recomposing or any_level_changed) {
-
-      for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-        (*f)->recompose_allocate (rl);
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_fill (state, rl, do_prolongate);
-        }
-        (*f)->recompose_free (rl);
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
-        }
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_sync (state, rl, do_prolongate);
-        }
-      } // for all grid functions of same vartype
-
-    } // if did_change
-  } // for all refinement levels
-}
-
-void dh::save_memory (bool do_prolongate)
-{
-  DECLARE_CCTK_PARAMETERS;
-
-  for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-    
-    (*f)->recompose_crop ();
-
-    bool any_level_changed = false;
-    for (int rl=0; rl<h.reflevels(); ++rl) {
-      // TODO: calculate this_level_changed only once
-      bool const this_level_changed = (*f)->recompose_did_change (rl);
-      any_level_changed |= this_level_changed;
-      if (not fast_recomposing or any_level_changed) {
-
-        (*f)->recompose_allocate (rl);
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_fill (state, rl, do_prolongate);
-        }
-        (*f)->recompose_free (rl);
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
-        }
-        for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_sync (state, rl, do_prolongate);
-        }
-
-      } // if did_change
-    } // for rl
-    
-  } // for gf
-}
 
 // Grid function management
 void dh::add (ggf* f)
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 322f7bd5a..eaa6bcbc8 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -118,8 +118,9 @@ public:                         // should be readonly
   int inner_buffer_width;       // buffer inside refined grids
   i2vect buffers;               // buffer outside refined grids
   
-  mboxes boxes;
-  mbases bases;
+  mboxes boxes;                 // grid hierarchy
+  mbases bases;                 // bounding boxes around the grid
+                                // hierarchy
   
   list<ggf*> gfs;               // list of all grid functions
   
@@ -137,7 +138,8 @@ public:
   int prolongation_stencil_size () const;
   
   // Modifiers
-  void recompose (const bool do_prolongate);
+  void regrid ();
+  void recompose (const int rl, const bool do_prolongate);
   
   // Grid function management
   void add (ggf* f);
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 74d718d2a..a78b35def 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -106,30 +106,6 @@ void ggf::recompose_crop ()
   } // for ml
 }
 
-bool ggf::recompose_did_change (const int rl) const
-{
-  // Find out whether this level changed
-  if ((int)storage.size() != h.mglevels()) return true;
-  for (int ml=0; ml<h.mglevels(); ++ml) {
-    if ((int)storage.at(ml).size() <= rl) return true;
-    if ((int)storage.at(ml).at(rl).size() != h.components(rl)) return true;
-    for (int c=0; c<h.components(rl); ++c) {
-      if ((int)storage.at(ml).at(rl).at(c).size() != timelevels(ml,rl)) {
-        return true;
-      }
-      for (int tl=0; tl<timelevels(ml,rl); ++tl) {
-        ibbox const & wantextent = d.boxes.at(ml).at(rl).at(c).exterior;
-        ibbox const & haveextent = storage.at(ml).at(rl).at(c).at(tl)->extent();
-        if (haveextent != wantextent) return true;
-        int const wantproc = h.proc(rl,c);
-        int const haveproc = storage.at(ml).at(rl).at(c).at(tl)->proc();
-        if (haveproc != wantproc) return true;
-      } // for tl
-    } // for c
-  } // for ml
-  return false;
-}
-
 void ggf::recompose_allocate (const int rl)
 {
   // Retain storage that might be needed
@@ -543,6 +519,14 @@ void ggf::ref_bnd_prolongate (comm_state& state,
   vector<int> tl2s;
   if (transport_operator != op_copy) {
     // Interpolation in time
+    if (not (timelevels(ml,rl) >= prolongation_order_time+1)) {
+      char * const fullname = CCTK_FullName (varindex);
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The variable \"%s\" has only %d active time levels, which is not enough for boundary prolongation of order %d",
+                  fullname ? fullname : "<unknown variable>",
+                  timelevels(ml,rl), prolongation_order_time);
+      free (fullname);
+    }
     assert (timelevels(ml,rl) >= prolongation_order_time+1);
     tl2s.resize(prolongation_order_time+1);
     for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 5a2bab361..ebc530ef9 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -91,7 +91,6 @@ public:
   void set_timelevels (int ml, int rl, int new_timelevels);
 
   void recompose_crop ();
-  bool recompose_did_change (int rl) const;
   void recompose_allocate (int rl);
   void recompose_fill (comm_state& state, int rl, bool do_prolongate);
   void recompose_free (int rl);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 6765e946e..24dd0d0e2 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -37,13 +37,17 @@ gh::gh (const vector<ivect> & reffacts_, const centering refcent_,
 gh::~gh () { }
 
 // Modifiers
-void gh::recompose (const mexts& exts,
-                    const rbnds& outer_bounds,
-                    const rprocs& procs,
-                    const bool do_prolongate)
+void gh::regrid (const mexts& exts,
+                 const rbnds& outer_bounds,
+                 const rprocs& procs)
 {
   DECLARE_CCTK_PARAMETERS;
   
+  // Save the old grid hierarchy
+  _oldextents = _extents;
+  _oldouter_boundaries = _outer_boundaries;
+  _oldprocessors = _processors;
+  
   _extents = exts;
   _outer_boundaries = outer_bounds;
   _processors = procs;
@@ -66,16 +70,64 @@ void gh::recompose (const mexts& exts,
   }
   
   // Recompose the other hierarchies
-  
   for (list<th*>::iterator t=ths.begin(); t!=ths.end(); ++t) {
-    (*t)->recompose();
+    (*t)->regrid();
   }
   
   for (list<dh*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
-    (*d)->recompose (do_prolongate);
+    (*d)->regrid();
+  }
+}
+
+void gh::recompose (const int rl,
+                    const bool do_prolongate)
+{
+  // Handle changes in number of mglevels
+  if (_oldextents.size() != _extents.size()) {
+    _oldextents.resize (_extents.size());
+  }
+  
+  if (level_did_change(rl)) {
+    
+    // Recompose the other hierarchies
+    for (list<dh*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
+      (*d)->recompose (rl, do_prolongate);
+    }
+    
+    // Overwrite old with new grid hierarchy
+    for (int ml=0; ml<mglevels(); ++ml) {
+      _oldextents.at(ml).resize (_extents.at(ml).size());
+      _oldextents.at(ml).at(rl) = _extents.at(ml).at(rl);
+    }
+    _oldouter_boundaries.resize (_outer_boundaries.size());
+    _oldouter_boundaries.at(rl) = _outer_boundaries.at(rl);
+    _oldprocessors.resize (_processors.size());
+    _oldprocessors.at(rl) = _processors.at(rl);
   }
 }
 
+bool gh::level_did_change (const int rl) const
+{
+  // Find out whether this level changed
+  if (_extents.size() != _oldextents.size()) return true;
+  for (int ml=0; ml<mglevels(); ++ml) {
+    assert (rl>=0 and rl<reflevels());
+    if (rl >= (int)_oldextents.at(ml).size()) return true;
+    if (_extents.at(ml).at(rl).size() != _oldextents.at(ml).at(rl).size()) {
+      return true;
+    }
+    for (int c=0; c<components(rl); ++c) {
+      if (_extents.at(ml).at(rl).at(c).size() !=
+          _oldextents.at(ml).at(rl).at(c).size())
+      {
+        return true;
+      }
+      if (_processors.at(rl).at(c) != _oldprocessors.at(rl).at(c)) return true;
+    } // for c
+  } // for ml
+  return false;
+}
+
 void gh::check_processor_number_consistency ()
 {
   for (int rl=0; rl<reflevels(); ++rl) {
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index f46fd5c34..4e343c533 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -62,6 +62,10 @@ private:
   rbnds _outer_boundaries;	// boundary descriptions of all grids
   rprocs _processors;		// processor numbers of all grids
 
+  mexts _oldextents;            // a copy, used during regridding
+  rbnds _oldouter_boundaries;
+  rprocs _oldprocessors;
+
   list<th*> ths;		// list of all time hierarchies
   list<dh*> dhs;		// list of all data hierarchies
   
@@ -76,10 +80,15 @@ public:
   virtual ~gh ();
   
   // Modifiers
-  void recompose (const mexts& exts, const rbnds& outer_bounds,
-		  const rprocs& procs,
+  void regrid (const mexts& exts,
+               const rbnds& outer_bounds,
+               const rprocs& procs);
+  void recompose (const int rl,
                   const bool do_prolongate);
-
+private:
+  bool level_did_change (const int rl) const;
+  
+public:
   const mexts & extents() const
   {
     return _extents;
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index 0751a03e1..bf69d1b5c 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -34,7 +34,7 @@ th::~th ()
 }
 
 // Modifiers
-void th::recompose ()
+void th::regrid ()
 {
   const int old_mglevels = times.size();
   times.resize(h.mglevels());
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 3c24c227c..9772afaa5 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -47,7 +47,7 @@ public:
   ~th ();
   
   // Modifiers
-  void recompose ();
+  void regrid ();
   
   // Time management
   CCTK_REAL get_time (const int rl, const int ml) const