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

8 files changed, 113 insertions, 112 deletions

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