CarpetLib: use map instead of set for dh::gfs

some routines assume that elements in gfs are sorted by variable index
when traversing the container (eg recompose_allocate does)

4 files changed, 53 insertions, 21 deletions

diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 02474b8b9..a888417b1 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -129,6 +129,19 @@ memoryof (set<T> const & c)
   return s;
 }
 
+template <class S, class T>
+size_t
+memoryof (map<S,T> const & c)
+{
+  size_t s = sizeof c;
+  for (typename map<S,T>::const_iterator i=c.begin(); i!=c.end(); ++i) {
+    // Assume that there are three pointers per list element, forming
+    // a tree structure
+    s += 3 * sizeof (void *) + memoryof(i->second);
+  }
+  return s;
+}
+
 template <class T>
 size_t
 memoryof (stack<T> const & c)
@@ -345,6 +358,7 @@ template size_t memoryof (set<dh*> const & l);
 template size_t memoryof (set<gh*> const & l);
 template size_t memoryof (set<gdata*> const & l);
 template size_t memoryof (set<ggf*> const & l);
+template size_t memoryof (map<int,ggf*> const & l);
 template size_t memoryof (set<th*> const & l);
 template size_t memoryof (stack<void*> const & s);
 template size_t memoryof (vector<bool> const & v);
@@ -419,6 +433,7 @@ template istream& input (istream& os, vector<vector<vector<region_t> > >& v);
 
 //template ostream& output (ostream& os, const list<ibbox>& l);
 //template ostream& output (ostream& os, const list<region_t>& l);
+template ostream& output (ostream& os, const pair<int const, ggf*>& p);
 #ifdef CARPET_ENABLE_BBOXSET2
 //template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,0> > >& m);
 //template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,1> > >& m);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 747761cd7..8f8c13345 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -326,6 +326,8 @@ template<class T> inline size_t memoryof (typename list<T>::const_iterator const
 
 template<class T> size_t memoryof (list<T> const & c) CCTK_ATTRIBUTE_PURE;
 template<class T> size_t memoryof (set<T> const & c) CCTK_ATTRIBUTE_PURE;
+template<class S, class T>
+  size_t memoryof (map<S,T> const & c) CCTK_ATTRIBUTE_PURE;
 template<class T> size_t memoryof (stack<T> const & c) CCTK_ATTRIBUTE_PURE;
 template<class T> size_t memoryof (vector<T> const & c) CCTK_ATTRIBUTE_PURE;
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index b66c80188..8e93b8705 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -2115,10 +2115,10 @@ regrid (bool const do_init)
     cout << "memoryof(dh.fast_boxes)=" << memoryof(fast_boxes) << eol;
     int gfcount = 0;
     size_t gfmemory = 0;
-    for (set<ggf*>::const_iterator gfi = gfs.begin(); gfi != gfs.end(); ++ gfi)
+    for (map<int,ggf*>::const_iterator gfi = gfs.begin(); gfi != gfs.end(); ++ gfi)
     {
       ++ gfcount;
-      gfmemory += memoryof(**gfi);
+      gfmemory += memoryof(*gfi->second);
     }
     cout << "#gfs=" << gfcount << eol;
     cout << "memoryof(gfs)=" << gfmemory << eol;
@@ -2200,51 +2200,51 @@ recompose (int const rl, bool const do_prolongate)
   static Timers::Timer timer ("CarpetLib::dh::recompose");
   timer.start ();
   
-  for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-    (*f)->recompose_crop ();
+  for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+    f->second->recompose_crop ();
   }
   
   if (combine_recompose) {
     // Recompose all grid functions of this refinement levels at once.
     // This may be faster, but requires more memory.
-    for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-      (*f)->recompose_allocate (rl);
+    for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      f->second->recompose_allocate (rl);
     }
     // TODO: If this works, rename do_prolongate to do_init here, and
     // remove the do_prolongate parameter from ggf::recompose_fill
 #if 0
     for (comm_state state; not state.done(); state.step()) {
-      for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-        (*f)->recompose_fill (state, rl, do_prolongate);
+      for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+        f->second->recompose_fill (state, rl, do_prolongate);
       }
     }
 #endif
     if (do_prolongate) {
       for (comm_state state; not state.done(); state.step()) {
-        for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-          (*f)->recompose_fill (state, rl, true);
+        for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+          f->second->recompose_fill (state, rl, true);
         }
       }
     }
-    for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-      (*f)->recompose_free_old (rl);
+    for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      f->second->recompose_free_old (rl);
     }
   } else {
     // Recompose the grid functions sequentially.  This may be slower,
     // but requires less memory.  This is the default.
-    for (set<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-      (*f)->recompose_allocate (rl);
+    for (map<int,ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      f->second->recompose_allocate (rl);
 #if 0
       for (comm_state state; not state.done(); state.step()) {
-        (*f)->recompose_fill (state, rl, do_prolongate);
+        f->second->recompose_fill (state, rl, do_prolongate);
       }
 #endif
       if (do_prolongate) {
         for (comm_state state; not state.done(); state.step()) {
-          (*f)->recompose_fill (state, rl, true);
+          f->second->recompose_fill (state, rl, true);
         }
       }
-      (*f)->recompose_free_old (rl);
+      f->second->recompose_free_old (rl);
     }
   }
   
@@ -2259,7 +2259,12 @@ dh::
 insert (ggf * const f)
 {
   CHECKPOINT;
-  gfs.insert (f);
+  assert (f);
+  assert (f->varindex >= 0);
+  assert (f->varindex < CCTK_NumVars());
+  const bool inserted =
+    gfs.insert (pair<int,ggf*> (f->varindex, f)).second;
+  assert (inserted);
 }
 
 void
@@ -2267,7 +2272,11 @@ dh::
 erase (ggf * const f)
 {
   CHECKPOINT;
-  gfs.erase (f);
+  assert (f);
+  assert (f->varindex >= 0);
+  assert (f->varindex < CCTK_NumVars());
+  const size_t erased = gfs.erase (f->varindex);
+  assert (erased == 1);
 }
 
 
@@ -2750,7 +2759,7 @@ output (ostream & os)
      << "gfs={";
   {
     bool isfirst = true;
-    for (set<ggf*>::const_iterator
+    for (map<int,ggf*>::const_iterator
            f = gfs.begin(); f != gfs.end(); ++ f, isfirst = false)
     {
       if (not isfirst) os << ",";
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index b558b013f..55b75218c 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -5,6 +5,7 @@
 #include <cstdlib>
 #include <iostream>
 #include <set>
+#include <map>
 #include <string>
 #include <vector>
 
@@ -248,7 +249,12 @@ public:                         // should be readonly
   level_mboxes level_boxes;     // grid hierarchy [ml][rl]
   fast_mboxes fast_boxes;       // grid hierarchy [ml][rl][p]
   
-  set<ggf*> gfs;                // all grid functions
+private:
+
+  // this needs to be sorted by varindex so that when iteration through the
+  // container in order with a forward iterator, vector leaders are processed
+  // first
+  map<int,ggf*> gfs;            // all grid functions
   
 public: