CarpetLib: Introduce pseudoregion class

A pseudoregion is almost a region; it is a bbox that lives on a
certain processor.  Pseudoregions are a compact way to store
information about what processors needs to send data to what other
processor during synchronisation or regridding.

Introduce a class pseudoregion.  Add additional pseudoregion fields to
the dh class.  Use these fields when communicating; this is faster
than walking the original data structures.

darcs-hash:20070204184339-dae7b-3bd90a48b77cf0623e3df3a7ca4edd9d5c7053f8.gz

4 files changed, 161 insertions, 11 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index fb3a98538..7ea31573e 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -542,6 +542,43 @@ void dh::trim_unsynced_boundaries (dh::dboxes & box,
   }
 }
 
+void
+dh::
+optimise_field (dboxes & box,
+                iblistvect const dboxes::* const field,
+                pvect dboxes::* const field_fast,
+                int const rl, int const c, int const ml)
+{
+  size_t num_regions = 0;
+  for (int cc=0; cc<h.components(rl); ++cc) {
+    num_regions += (box.*field).AT(cc).size();
+  }
+  assert ((box.*field_fast).empty());
+  (box.*field_fast).reserve (num_regions);
+  for (int cc=0; cc<h.components(rl); ++cc) {
+    for (iblist::const_iterator
+           r=(box.*field).AT(cc).begin(); r!=(box.*field).AT(cc).end(); ++r)
+    {
+      pseudoregion pr;
+      pr.extent = * r;
+      pr.processor = cc;
+      (box.*field_fast).push_back (pr);
+    }
+  }
+  assert ((box.*field_fast).size() == num_regions);
+}
+
+void
+dh::
+optimise_fields (dboxes & box,
+                 int const rl, int const c, int const ml)
+{
+  optimise_field (box, &dboxes::recv_ref_fine, &dboxes::recv_ref_fine_fast, rl, c, ml);
+  optimise_field (box, &dboxes::recv_ref_coarse, &dboxes::recv_ref_coarse_fast, rl, c, ml);
+  optimise_field (box, &dboxes::recv_sync, &dboxes::recv_sync_fast, rl, c, ml);
+  optimise_field (box, &dboxes::recv_ref_bnd_coarse, &dboxes::recv_ref_bnd_coarse_fast, rl, c, ml);
+}
+
 void dh::check_bboxes (dh::dboxes & box,
                        int const rl, int const c, int const ml)
 {
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index bddc4360b..86e3c4f1a 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -17,6 +17,17 @@ using namespace std;
 
 
 
+// A pseudoregion is almost a region; it is a bbox that lives on a
+// certain processor.  Pseudoregions are a compact way to store
+// information about what processors needs to send data to what other
+// processor during synchronisation or regridding.
+struct pseudoregion {
+  ibbox extent;
+  int processor;
+};
+
+
+
 // Forward declaration
 class ggf;
 class dh;
@@ -34,6 +45,9 @@ public:
   typedef list<ibbox>    iblist;
   typedef vector<iblist> iblistvect; // vector of lists
 
+  typedef vector <pseudoregion> pvect;
+
+
 
   // in here, the term "boundary" means both ghost zones and
   // refinement boundaries, but does not refer to outer (physical)
@@ -57,13 +71,17 @@ public:
     iblistvect send_ref_fine;
     iblistvect send_ref_coarse;
     iblistvect recv_ref_fine;
+    pvect recv_ref_fine_fast;
     iblistvect recv_ref_coarse;
+    pvect recv_ref_coarse_fast;
     iblistvect send_sync;       // send while syncing
     iblistvect send_ref_bnd_fine; // sent to finer grids
     
     ibset boundaries;           // boundaries
     iblistvect recv_sync;       // received while syncing
+    pvect recv_sync_fast;
     iblistvect recv_ref_bnd_coarse; // received from coarser grids
+    pvect recv_ref_bnd_coarse_fast;
     ibset sync_not;             // not received while syncing (outer
                                 // boundary of that level)
     ibset recv_not;             // not received while syncing or
@@ -100,6 +118,12 @@ private:
   void setup_refinement_prolongation_boxes (dboxes & b, int rl, int c, int ml);
   void setup_refinement_boundary_prolongation_boxes (dboxes & b, int rl, int c, int ml);
   void setup_refinement_restriction_boxes (dboxes & b, int rl, int c, int ml);
+  void optimise_field (dboxes & b,
+                       iblistvect const dboxes::* field,
+                       pvect dboxes::* field_fast,
+                       int rl, int c, int ml);
+  void optimise_fields (dboxes & b,
+                        int rl, int c, int ml);
   void trim_unsynced_boundaries (dboxes & b, int rl, int c, int ml);
   void do_output_bboxes (dboxes & b, int rl, int c, int ml);
   void check_bboxes (dboxes & b, int rl, int c, int ml);
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 0778d83fe..7009ed0f9 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -375,16 +375,51 @@ void ggf::copycat (comm_state& state,
   // walk all components
   static Timer copycat1 ("copycat_listvect_1");
   copycat1.start ();
-  for (int c2=0; c2<h.components(rl2); ++c2) {
-    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+  iblistvect::const_iterator recvi =
+    (d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_listvect).begin();
+  gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
+  cdata::const_iterator srci = storage.AT(ml2).AT(rl2).begin();
+  int const nc = h.components(rl2);
+  for (int c2=0; c2<nc; ++c2) {
+    iblist const & recv = (*recvi);
+    gdata const * const src = (*srci).AT(tl2);
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(); r!=recv.end(); ++r) {
       // (use the send boxes for communication)
       // copy the content
-      gdata* const dst = storage.at(ml1).at(rl1).at(c1).at(tl1);
-      gdata* const src = storage.at(ml2).at(rl2).at(c2).at(tl2);
       dst->copy_from(state, src, *r);
     }
+    ++ recvi;
+    ++ srci;
+  }
+  copycat1.stop (0);
+}
+
+// Copy regions
+void ggf::copycat (comm_state& state,
+                   int tl1, int rl1, int c1, int ml1,
+                   const pvect dh::dboxes::* recv_pvect,
+                   int tl2, int rl2, int ml2)
+{
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl1>=0 and tl1<timelevels(ml1,rl1));
+  assert (           ml2<h.mglevels());
+  assert (rl2>=0 and rl2<h.reflevels());
+  assert (tl2>=0 and tl2<timelevels(ml2,rl2));
+  // walk all components
+  static Timer copycat1 ("copycat_pvect_1");
+  copycat1.start ();
+  gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
+  cdata const & srcs = storage.AT(ml2).AT(rl2);
+  pvect const & prs = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_pvect;
+  for (pvect::const_iterator ipr=prs.begin(); ipr!=prs.end(); ++ipr) {
+    pseudoregion const & pr = * ipr;
+    ibbox const & r = pr.extent;
+    int const c2 = pr.processor;
+    gdata const * const src = srcs.AT(c2).AT(tl2);
+    dst->copy_from(state, src, r);
   }
   copycat1.stop (0);
 }
@@ -509,6 +544,45 @@ void ggf::intercat (comm_state& state,
   intercat1.stop (0);
 }
 
+// Interpolate regions
+void ggf::intercat (comm_state& state,
+                    int tl1, int rl1, int c1, int ml1,
+                    const pvect dh::dboxes::* recv_pvect,
+                    const vector<int> tl2s, int rl2, int ml2,
+                    const CCTK_REAL time)
+{
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl1>=0 and tl1<timelevels(ml1,rl1));
+  assert (rl2>=0 and rl2<h.reflevels());
+  assert (ml2>=0 and ml2<h.mglevels());
+  vector<CCTK_REAL> times(tl2s.size());
+  for (int i=0; i<(int)tl2s.size(); ++i) {
+    assert (tl2s.AT(i)>=0 and tl2s.AT(i)<timelevels(ml2,rl2));
+    times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
+  }
+  // walk all components
+  static Timer intercat1 ("intercat_pvect_1");
+  intercat1.start ();
+  int const pos = d.prolongation_order_space;
+  int const pot = transport_operator == op_copy ? 0 : prolongation_order_time;
+  gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
+  cdata const & srcs = storage.AT(ml2).AT(rl2);
+  vector<const gdata*> gsrcs(tl2s.size());
+  pvect const & prs = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_pvect;
+  for (pvect::const_iterator ipr=prs.begin(); ipr!=prs.end(); ++ipr) {
+    pseudoregion const & pr = * ipr;
+    ibbox const & r = pr.extent;
+    int const c2 = pr.processor;
+    for (int i=0; i<(int)gsrcs.size(); ++i) {
+      gsrcs.AT(i) = srcs.AT(c2).AT(tl2s.AT(i));
+    }
+    dst->interpolate_from (state, gsrcs, times, r, time, pos, pot);
+  }
+  intercat1.stop (0);
+}
+
 
 
 // Copy a component from the next time level
@@ -530,7 +604,7 @@ void ggf::sync (comm_state& state, int tl, int rl, int c, int ml)
   static Timer timer ("sync");
   timer.start ();
   copycat (state,
-           tl,rl,c,ml, &dh::dboxes::recv_sync,
+           tl,rl,c,ml, &dh::dboxes::recv_sync_fast,
       	   tl,rl,  ml);
   timer.stop (0);
 }
@@ -565,7 +639,7 @@ void ggf::ref_bnd_prolongate (comm_state& state,
     tl2s.AT(0) = 0;
   }
   intercat (state,
-            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse,
+            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse_fast,
             tl2s,rl-1,  ml, time);
   timer.stop (0);
 }
@@ -613,11 +687,11 @@ void ggf::ref_restrict (comm_state& state,
   assert (abs(t.get_time(rl,ml) - t.get_time(rl+1,ml))
 	  <= 1.0e-8 * abs(t.get_time(rl,ml)));
   if (transport_operator == op_none) return;
-  const vector<int> tl2s(1,tl);
   static Timer timer ("ref_restrict");
   timer.start ();
+  vector<int> const tl2s(1,tl);
   intercat (state,
-            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_fine,
+            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_fine_fast,
 	    tl2s,rl+1,  ml, time);
   timer.stop (0);
 }
@@ -629,15 +703,15 @@ void ggf::ref_prolongate (comm_state& state,
 {
   assert (rl>=1);
   if (transport_operator == op_none) return;
-  vector<int> tl2s;
   static Timer timer ("ref_prolongate");
   timer.start ();
+  vector<int> tl2s;
   // Interpolation in time
   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;
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s.AT(i) = i;
   intercat (state,
-            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_coarse,
+            tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_coarse_fast,
 	    tl2s,rl-1,  ml, time);
   timer.stop (0);
 }
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 30752a422..3e113ee01 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -33,6 +33,8 @@ class ggf {
   typedef list<ibbox>    iblist;
   typedef vector<iblist> iblistvect;
   
+  typedef vector <pseudoregion> pvect;
+  
   typedef gdata*        tdata;  // data ...
   typedef vector<tdata> fdata;  // ... for each time level
   typedef vector<fdata> cdata;  // ... for each component
@@ -134,6 +136,12 @@ protected:
                 int tl1, int rl1, int c1, int ml1,
 		const iblistvect dh::dboxes::* recv_listvect,
 		int tl2, int rl2, int ml2);
+
+  // Copy regions
+  void copycat (comm_state& state,
+                int tl1, int rl1, int c1, int ml1,
+		const pvect dh::dboxes::* recv_pvect,
+                int tl2, int rl2, int ml2);
   
   // Interpolate a region
   void intercat (comm_state& state,
@@ -156,6 +164,13 @@ protected:
 		 const vector<int> tl2s, int rl2, int ml2,
 		 CCTK_REAL time);
 
+  // Interpolate regions
+  void intercat (comm_state& state,
+                 int tl1, int rl1, int c1, int ml1,
+		 const pvect dh::dboxes::* recv_pvect,
+		 const vector<int> tl2s, int rl2, int ml2,
+		 CCTK_REAL time);
+
 
 
 public: