Remove the parameters Carpet::prolongate_initial_data; this is now

Remove the parameters Carpet::prolongate_initial_data; this is now
always done.

Remove arguments initialise_from and do_prolongate from Regrid().

Regridding is now done in level mode instead of meta mode.
Furthermore, CarpetRegrid is called in singlemape mode.

darcs-hash:20040418112943-07bb3-2e392df1737ab75f3f0d553bb53bde2ed41f8773.gz

7 files changed, 89 insertions, 71 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 1cd9772e7..f7b7b88e5 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.50 2004/04/07 16:56:42 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.51 2004/04/18 13:29:43 schnetter Exp $
 
 #include <assert.h>
 
@@ -31,7 +31,7 @@ dh<D>::dh (gh<D>& h,
   assert (buffer_width>=0);
   h.add(this);
   CHECKPOINT;
-  recompose (0, true);
+  recompose ();
 }
 
 // Destructors
@@ -51,7 +51,7 @@ int dh<D>::prolongation_stencil_size () const {
 
 // Modifiers
 template<int D>
-void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
+void dh<D>::recompose () {
   DECLARE_CCTK_PARAMETERS;
   
   CHECKPOINT;
@@ -574,7 +574,7 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
   
   for (typename list<ggf<D>*>::iterator f=gfs.begin();
        f!=gfs.end(); ++f) {
-    (*f)->recompose (initialise_from, do_prolongate);
+    (*f)->recompose ();
   }
 }
 
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index d525a0be4..c78d31291 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.17 2004/03/03 16:20:19 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.18 2004/04/18 13:29:43 schnetter Exp $
 
 #ifndef DH_HH
 #define DH_HH
@@ -115,7 +115,7 @@ public:
   int prolongation_stencil_size () const;
   
   // Modifiers
-  void recompose (const int initialise_from, const bool do_prolongate);
+  void recompose ();
   
   // Grid function management
   void add (ggf<D>* f);
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 462c179b8..36b601850 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.17 2004/03/23 19:30:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.18 2004/04/18 13:29:43 schnetter Exp $
 
 #include <assert.h>
 
@@ -23,7 +23,7 @@ gf<T,D>::gf (const int varindex, const operator_type transport_operator,
            t, d, tmin, tmax, prolongation_order_time,
            vectorlength, vectorindex, vectorleader)
 {
-  this->recompose (0, true);
+  this->recompose ();
 }
 
 // Destructors
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 38657d367..6b78c43dd 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.36 2004/03/23 19:30:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.37 2004/04/18 13:28:25 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -55,7 +55,7 @@ bool ggf<D>::operator== (const ggf<D>& f) const {
 
 // Modifiers
 template<int D>
-void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
+void ggf<D>::recompose () {
   
   // TODO: restructure storage only when needed
   
@@ -66,7 +66,7 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
   storage.resize(tmax-tmin+1);
   for (int tl=tmin; tl<=tmax; ++tl) {
     storage.at(tl-tmin).resize(h.reflevels());
-    for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
       storage.at(tl-tmin).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));
@@ -78,53 +78,73 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
   } // for tl
   
   // Initialise the new storage
-  for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
-    for (int tl=tmin; tl<=tmax; ++tl) {
-      for (int c=0; c<h.components(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);
-	  
-      	  // Allocate storage
-      	  storage.at(tl-tmin).at(rl).at(c).at(ml)->allocate
-      	    (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c));
-	  
-          if (do_prolongate) {
-            // Initialise from coarser level, if possible
-            // TODO: init only un-copied regions
-            if (rl>0) {
-              for (comm_state<D> state; !state.done(); state.step()) {
-                const CCTK_REAL time = t.time(tl,rl,ml);
-                ref_prolongate (state,tl,rl,c,ml,time);
-              }
-            } // if rl
-          } // if do_prolongate
-          
-          // Copy from old storage, if possible
-          // todo: copy only from interior regions?
-          if (rl<(int)oldstorage.at(tl-tmin).size()) {
-            for (comm_state<D> state; !state.done(); state.step()) {
+  // TODO: overlap this initialisation for all variables
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    bool firsttime=true;
+    for (comm_state<D> state; !state.done(); state.step(), firsttime=false) {
+      for (int tl=tmin; tl<=tmax; ++tl) {
+        for (int c=0; c<h.components(rl); ++c) {
+          for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
+            
+            if (firsttime) {
+              storage.at(tl-tmin).at(rl).at(c).at(ml) = typed_data(tl,rl,c,ml);
+              
+              // Allocate storage
+              storage.at(tl-tmin).at(rl).at(c).at(ml)->allocate
+                (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c));
+            }
+            
+            // Find out which regions need to be 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);
+            
+            // 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()) {
-                  const ibbox ovlp =
-                    (d.boxes.at(rl).at(c).at(ml).exterior
-                     & oldstorage.at(tl-tmin).at(rl).at(cc).at(ml)->extent());
-                  storage.at(tl-tmin).at(rl).at(c).at(ml)->copy_from
-                    (state, oldstorage.at(tl-tmin).at(rl).at(cc).at(ml), ovlp);
+                  // 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 (typename 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
               } // for cc
-            } // for step
-          } // if rl
-          
-	} // for ml
-      } // for c
-      
-    } // for tl
+            } // if rl
+            work.normalize();
+            
+            // Initialise from coarser level, if possible
+            if (rl>0) {
+              if (transport_operator != op_none) {
+                assert (tmax-tmin+1 >= prolongation_order_time+1);
+                vector<int> tls(prolongation_order_time+1);
+                for (int i=0; i<=prolongation_order_time; ++i) {
+                  tls.at(i) = tmax - i;
+                }
+                vector<const gdata<D>*> gsrcs(tls.size());
+                vector<CCTK_REAL> times(tls.size());
+                for (int i=0; i<(int)gsrcs.size(); ++i) {
+                  gsrcs.at(i) = storage.at(tls.at(i)-tmin).at(rl-1).at(c).at(ml);
+                  times.at(i) = t.time(tls.at(i),rl-1,ml);
+                }
+                const CCTK_REAL time = t.time(tl,rl,ml);
+                for (typename ibset::const_iterator r=work.begin(); r!=work.end(); ++r) {
+                  storage.at(tl-tmin).at(rl).at(c).at(ml)->interpolate_from (state, gsrcs, times, *r, time, d.prolongation_order_space, prolongation_order_time);
+                }
+              } // if transport_operator
+            } // if rl
+            
+          } // for ml
+        } // for c
+      } // for tl
+    } // for step
   } // for rl
   
   // Delete old storage
   for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int rl=initialise_from; rl<(int)oldstorage.at(tl-tmin).size(); ++rl) {
+    for (int rl=0; rl<(int)oldstorage.at(tl-tmin).size(); ++rl) {
       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);
@@ -134,21 +154,19 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
   } // for tl
   
   for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
       
       // Set boundaries
       for (int c=0; c<h.components(rl); ++c) {
       	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
           
-          if (do_prolongate) {
-            // TODO: assert that reflevel 0 boundaries are copied
-            if (rl>0) {
-              for (comm_state<D> state; !state.done(); state.step()) {
-                const CCTK_REAL time = t.time(tl,rl,ml);
-                ref_bnd_prolongate (state,tl,rl,c,ml,time);
-              }
-            } // if rl
-          }
+          // TODO: assert that reflevel 0 boundaries are copied
+          if (rl>0) {
+            for (comm_state<D> state; !state.done(); state.step()) {
+              const CCTK_REAL time = t.time(tl,rl,ml);
+              ref_bnd_prolongate (state,tl,rl,c,ml,time);
+            }
+          } // if rl
           
           for (comm_state<D> state; !state.done(); state.step()) {
             sync (state,tl,rl,c,ml);
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 64998c8a2..d228a863e 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.20 2004/03/23 12:40:27 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.21 2004/04/18 13:29:43 schnetter Exp $
 
 #ifndef GGF_HH
 #define GGF_HH
@@ -87,7 +87,7 @@ public:
 
 
   // Modifiers
-  void recompose (const int initialise_from, const bool do_prolongate);
+  void recompose ();
 
   // Cycle the time levels by rotating the data sets
   void cycle (int rl, int c, int ml);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 172da013d..b68be34d8 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.26 2004/04/07 16:59:47 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.28 2004/04/18 13:29:43 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -36,9 +36,7 @@ gh<D>::~gh () { }
 template<int D>
 void gh<D>::recompose (const rexts& exts,
                        const rbnds& outer_bounds,
-		       const rprocs& procs,
-                       const int initialise_from,
-                       const bool do_prolongate)
+		       const rprocs& procs)
 {
   DECLARE_CCTK_PARAMETERS;
   
@@ -96,7 +94,11 @@ void gh<D>::recompose (const rexts& exts,
   // Check base grid extent
   if (reflevels()>0) {
     for (int c=0; c<components(0); ++c) {
+      // TODO: put the check back in, taking outer boundaries into
+      // account
+#if 0
       assert (extents.at(0).at(c).at(0).is_contained_in(baseextent));
+#endif
     }
   }
   
@@ -169,7 +171,7 @@ void gh<D>::recompose (const rexts& exts,
   }
   
   for (typename list<dh<D>*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
-    (*d)->recompose (initialise_from, do_prolongate);
+    (*d)->recompose ();
   }
 }
 
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index babc85cfb..69f9f31c0 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.16 2004/03/23 19:30:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.17 2004/04/18 13:29:43 schnetter Exp $
 
 #ifndef GH_HH
 #define GH_HH
@@ -85,9 +85,7 @@ public:
   
   // Modifiers
   void recompose (const rexts& exts, const rbnds& outer_bounds,
-		  const rprocs& procs,
-                  const int initialise_from,
-                  const bool do_prolongate);
+		  const rprocs& procs);
   
   // Accessors
   int reflevels () const {