Do not initialise the finer levels of the grid functions at initial

Do not initialise the finer levels of the grid functions at initial
time through prolongation.

darcs-hash:20030502122312-07bb3-ad2023a18bb9177cae3fbbe0df5aa4401f0a0259.gz

6 files changed, 55 insertions, 41 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 60a7500d1..61138d17e 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.30 2003/04/30 12:39:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.31 2003/05/02 14:23:12 schnetter Exp $
 
 #include <assert.h>
 
@@ -27,7 +27,7 @@ dh<D>::dh (gh<D>& h, const ivect& lghosts, const ivect& ughosts,
   assert (all(lghosts>=0 && ughosts>=0));
   h.add(this);
   CHECKPOINT;
-  recompose();
+  recompose(false);
 }
 
 // Destructors
@@ -47,7 +47,7 @@ int dh<D>::prolongation_stencil_size () const {
 
 // Modifiers
 template<int D>
-void dh<D>::recompose () {
+void dh<D>::recompose (const int initialise_upto) {
   DECLARE_CCTK_PARAMETERS;
   
   CHECKPOINT;
@@ -478,7 +478,7 @@ void dh<D>::recompose () {
   
   for (typename list<ggf<D>*>::iterator f=gfs.begin();
        f!=gfs.end(); ++f) {
-    (*f)->recompose();
+    (*f)->recompose(initialise_upto);
   }
 }
 
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 115ead405..98c7c75e3 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.13 2003/03/26 17:34:43 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.14 2003/05/02 14:23:12 schnetter Exp $
 
 #ifndef DH_HH
 #define DH_HH
@@ -114,7 +114,7 @@ public:
   int prolongation_stencil_size () const;
   
   // Modifiers
-  void recompose ();
+  void recompose (const int initialise_upto);
   
   // Grid function management
   void add (ggf<D>* f);
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 7ab793e0f..067542fc8 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.23 2003/03/28 10:11:54 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.24 2003/05/02 14:23:12 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -50,7 +50,7 @@ bool ggf<D>::operator== (const ggf<D>& f) const {
 
 // Modifiers
 template<int D>
-void ggf<D>::recompose () {
+void ggf<D>::recompose (const int initialise_upto) {
   
   // Retain storage that might be needed
   fdata oldstorage = storage;
@@ -82,26 +82,31 @@ void ggf<D>::recompose () {
       	  storage[tl-tmin][rl][c][ml]->allocate
       	    (d.boxes[rl][c][ml].exterior, h.proc(rl,c));
 	  
-      	  // Initialise from coarser level, if possible
-          // TODO: init only un-copied regions
-      	  if (rl>0) {
-	    const CCTK_REAL time = t.time(tl,rl,ml);
-      	    ref_prolongate (tl,rl,c,ml,time);
-      	  } // if rl
-	  
-      	  // Copy from old storage, if possible
-      	  if (rl<(int)oldstorage[tl-tmin].size()) {
-      	    for (int cc=0; cc<(int)oldstorage[tl-tmin][rl].size(); ++cc) {
-      	      if (ml<(int)oldstorage[tl-tmin][rl][cc].size()) {
-      		const ibbox ovlp =
-      		  (d.boxes[rl][c][ml].exterior
-      		   & oldstorage[tl-tmin][rl][cc][ml]->extent());
-      		storage[tl-tmin][rl][c][ml]->copy_from
-      		  (oldstorage[tl-tmin][rl][cc][ml], ovlp);
-	      } // if ml
-	    } // for cc
-	  } // if rl
-	  
+          // Initialise only if desired
+          if (initialise_upto >= 0 && rl <= initialise_upto) {
+            
+            // Initialise from coarser level, if possible
+            // TODO: init only un-copied regions
+            if (rl>0) {
+              const CCTK_REAL time = t.time(tl,rl,ml);
+              ref_prolongate (tl,rl,c,ml,time);
+            } // if rl
+            
+            // Copy from old storage, if possible
+            if (rl<(int)oldstorage[tl-tmin].size()) {
+              for (int cc=0; cc<(int)oldstorage[tl-tmin][rl].size(); ++cc) {
+                if (ml<(int)oldstorage[tl-tmin][rl][cc].size()) {
+                  const ibbox ovlp =
+                    (d.boxes[rl][c][ml].exterior
+                     & oldstorage[tl-tmin][rl][cc][ml]->extent());
+                  storage[tl-tmin][rl][c][ml]->copy_from
+                    (oldstorage[tl-tmin][rl][cc][ml], ovlp);
+                } // if ml
+              } // for cc
+            } // if rl
+            
+          } // if initialise
+          
 	} // for ml
       } // for c
       
@@ -125,12 +130,19 @@ void ggf<D>::recompose () {
       // Set boundaries
       for (int c=0; c<h.components(rl); ++c) {
       	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	  sync (tl,rl,c,ml);
-      	  // TODO: assert that reflevel 0 boundaries are copied
-      	  if (rl>0) {
-	    const CCTK_REAL time = t.time(tl,rl,ml);
-      	    ref_bnd_prolongate (tl,rl,c,ml,time);
-      	  } // if rl
+          
+          // Initialise only if desired
+          if (initialise_upto >= 0 && rl <= initialise_upto) {
+            
+            sync (tl,rl,c,ml);
+            // TODO: assert that reflevel 0 boundaries are copied
+            if (rl>0) {
+              const CCTK_REAL time = t.time(tl,rl,ml);
+              ref_bnd_prolongate (tl,rl,c,ml,time);
+            } // if rl
+            
+          } // if initialise
+            
       	} // for ml
       } // for c
       
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 6c467632e..7b58c9117 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.13 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.14 2003/05/02 14:23:12 schnetter Exp $
 
 #ifndef GGF_HH
 #define GGF_HH
@@ -79,7 +79,7 @@ public:
 
 
   // Modifiers
-  void recompose ();
+  void recompose (const int initialise_upto = -1);
 
   // 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 c34abdb35..00cbad2d7 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.19 2003/04/30 12:39:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.20 2003/05/02 14:23:12 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -36,7 +36,8 @@ gh<D>::~gh () { }
 template<int D>
 void gh<D>::recompose (const rexts& exts,
                        const rbnds& outer_bounds,
-		       const rprocs& procs)
+		       const rprocs& procs,
+                       const int initialise_upto)
 {
   DECLARE_CCTK_PARAMETERS;
   
@@ -160,7 +161,7 @@ void gh<D>::recompose (const rexts& exts,
   }
   
   for (typename list<dh<D>*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
-    (*d)->recompose();
+    (*d)->recompose(initialise_upto);
   }
 }
 
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 3a2fb649c..c748e93ee 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.12 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.13 2003/05/02 14:23:12 schnetter Exp $
 
 #ifndef GH_HH
 #define GH_HH
@@ -85,7 +85,8 @@ public:
   
   // Modifiers
   void recompose (const rexts& exts, const rbnds& outer_bounds,
-		  const rprocs& procs);
+		  const rprocs& procs,
+                  const int initialise_upto = -1);
   
   // Helpers
   cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts,