Further changes to make Carpet work with multiple multigrid levels.

Further changes to make Carpet work with multiple multigrid levels.
Fixed one other nasty bug in a prolongation operator.

Optimised last remaining non-optimised prolongation operator.

darcs-hash:20020109164539-07bb3-6ebee1b591a732eb826557128a2a0bce38151ed1.gz

16 files changed, 291 insertions, 197 deletions

diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index 50b57ce38..5fc0b056b 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -9,7 +9,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.6 2001/12/09 16:41:52 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.7 2002/01/09 17:45:39 schnetter Exp $";
 
 
 
@@ -43,7 +43,6 @@ namespace Carpet {
     
     assert (group<(int)arrdata.size());
     for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
-//       if (num_tl>1 && reflevel>0) {
       if (reflevel>0) {
 	for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
 	  arrdata[group].data[var]->ref_bnd_prolongate
diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc
index 539c48bec..6c68cbd82 100644
--- a/Carpet/Carpet/src/Evolve.cc
+++ b/Carpet/Carpet/src/Evolve.cc
@@ -9,7 +9,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.6 2001/12/14 16:39:06 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.7 2002/01/09 17:45:39 schnetter Exp $";
 
 
 
@@ -66,76 +66,92 @@ namespace Carpet {
       
       Waypoint ("Evolving iteration %d...", cgh->cctk_iteration);
       
-      BEGIN_REFLEVEL_LOOP(cgh) {
-	if ((cgh->cctk_iteration-1) % (maxreflevelfact/reflevelfact) == 0) {
+      BEGIN_MGLEVEL_LOOP(cgh) {
+	if ((cgh->cctk_iteration-1) % mglevelfact == 0) {
 	  
-	  // Cycle time levels
-	  CycleTimeLevels (cgh);
+	  Waypoint ("Evolving multigrid level %d...", mglevel);
 	  
-	  // Advance level times
-	  tt->advance_time (reflevel, mglevel);
-	  tt0->advance_time (reflevel, mglevel);
-	  for (int group=0; group<CCTK_NumGroups(); ++group) {
-	    switch (CCTK_GroupTypeI(group)) {
-	    case CCTK_SCALAR:
-	      break;
-	    case CCTK_ARRAY:
-	      arrdata[group].tt->advance_time (reflevel, mglevel);
-	      break;
-	    case CCTK_GF:
-	      break;
-	    default:
-	      abort();
+	  BEGIN_REFLEVEL_LOOP(cgh) {
+	    const int do_every = mglevelfact * maxreflevelfact/reflevelfact;
+	    if ((cgh->cctk_iteration-1) % do_every == 0) {
+	      
+	      // Cycle time levels
+	      CycleTimeLevels (cgh);
+	      
+	      // Advance level times
+	      tt->advance_time (reflevel, mglevel);
+	      tt0->advance_time (reflevel, mglevel);
+	      for (int group=0; group<CCTK_NumGroups(); ++group) {
+		switch (CCTK_GroupTypeI(group)) {
+		case CCTK_SCALAR:
+		  break;
+		case CCTK_ARRAY:
+		  arrdata[group].tt->advance_time (reflevel, mglevel);
+		  break;
+		case CCTK_GF:
+		  break;
+		default:
+		  abort();
+		}
+	      }
+	      
+	      // Checking
+	      CalculateChecksums (cgh, allbutcurrenttime);
+	      Poison (cgh, currenttimebutnotifonly);
+	      
+	      // Evolve
+	      Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
+	      CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
+	      Waypoint ("%*sScheduling EVOL", 2*reflevel, "");
+	      CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
+	      Waypoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
+	      CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
+	      
+	      // Checking
+	      PoisonCheck (cgh, currenttimebutnotifonly);
+	      
+	      // Recompose grid hierarchy
+	      Recompose (cgh);
+	      
 	    }
-	  }
-	  
-	  // Checking
-	  CalculateChecksums (cgh, allbutcurrenttime);
-	  Poison (cgh, currenttimebutnotifonly);
-	  
-	  // Evolve
-	  Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
-	  CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
-	  Waypoint ("%*sScheduling EVOL", 2*reflevel, "");
-	  CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
-	  Waypoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
-	  CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
-	  
-	  // Checking
-	  PoisonCheck (cgh, currenttimebutnotifonly);
-	  
-	  // Recompose grid hierarchy
-	  Recompose (cgh);
+	  } END_REFLEVEL_LOOP(cgh);
 	  
 	}
-      } END_REFLEVEL_LOOP(cgh);
+      } END_MGLEVEL_LOOP(cgh);
       
-      BEGIN_REVERSE_REFLEVEL_LOOP(cgh) {
-	if (cgh->cctk_iteration % (maxreflevelfact/reflevelfact) == 0) {
-	  
-	  // Restrict
-	  Restrict (cgh);
+      BEGIN_MGLEVEL_LOOP(cgh) {
+	if (cgh->cctk_iteration % mglevelfact == 0) {
 	  
-	  // Checking
-	  CalculateChecksums (cgh, currenttime);
-	  
-	  // Waypoint
-	  Waypoint ("%*sScheduling CHECKPOINT", 2*reflevel, "");
-	  CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction);
-	  
-	  // Analysis
-	  Waypoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
-	  CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
-	  
-	  // Output
-	  Waypoint ("%*sOutputGH", 2*reflevel, "");
-	  CCTK_OutputGH (cgh);
-	  
-	  // Checking
-	  CheckChecksums (cgh, alltimes);
+	  BEGIN_REVERSE_REFLEVEL_LOOP(cgh) {
+	    const int do_every = mglevelfact * maxreflevelfact/reflevelfact;
+	    if (cgh->cctk_iteration % do_every == 0) {
+	      
+	      // Restrict
+	      Restrict (cgh);
+	      
+	      // Checking
+	      CalculateChecksums (cgh, currenttime);
+	      
+	      // Waypoint
+	      Waypoint ("%*sScheduling CHECKPOINT", 2*reflevel, "");
+	      CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction);
+	      
+	      // Analysis
+	      Waypoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
+	      CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
+	      
+	      // Output
+	      Waypoint ("%*sOutputGH", 2*reflevel, "");
+	      CCTK_OutputGH (cgh);
+	      
+	      // Checking
+	      CheckChecksums (cgh, alltimes);
+	      
+	    }
+	  } END_REVERSE_REFLEVEL_LOOP(cgh);
 	  
 	}
-      } END_REVERSE_REFLEVEL_LOOP(cgh);
+      } END_MGLEVEL_LOOP(cgh);
       
     } // main loop
     
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 76417d954..a10a3bb3d 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -9,7 +9,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.3 2001/12/18 10:29:35 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.4 2002/01/09 17:45:39 schnetter Exp $";
 
 
 
@@ -46,6 +46,8 @@ namespace Carpet {
     CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction);
     CCTKi_FinaliseParamWarn();
     
+    Waypoint ("Initialising iteration %d...", cgh->cctk_iteration);
+    
     BEGIN_REFLEVEL_LOOP(cgh) {
       
       // Checking
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 3792ee082..544e2406e 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -15,7 +15,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.15 2002/01/09 13:56:25 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.16 2002/01/09 17:45:39 schnetter Exp $";
 
 
 
@@ -83,6 +83,8 @@ namespace Carpet {
   void RegisterRecomposeRegions (const gh<dim>::rexts& bbsss,
 				 const gh<dim>::rprocs& pss)
   {
+    assert (mglevel == 0);
+    
     // Check the regions
     CheckRegions (bbsss, pss);
     // Save the region information for the next regridding
@@ -96,6 +98,9 @@ namespace Carpet {
   void Recompose (const cGH* cgh)
   {
     assert (component == -1);
+    
+    if (mglevel > 0) return;
+    
     Waypoint ("%*sRecompose", 2*reflevel, "");
     
     // Check whether to recompose
@@ -212,7 +217,9 @@ namespace Carpet {
     if (verbose) {
       cout << endl;
       cout << "New bounding boxes";
-      if (descr) {
+      if (!descr) {
+	cout << " for grid functions";
+      } else {
 	if (strlen(descr)) {
 	  cout << " for group " << descr;
 	} else {
@@ -230,7 +237,9 @@ namespace Carpet {
       }
       cout << endl;
       cout << "New processor distribution";
-      if (descr) {
+      if (!descr) {
+	cout << " for grid functions";
+      } else {
 	if (strlen(descr)) {
 	  cout << " for group " << descr;
 	} else {
diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc
index ad264a346..7dfa8b372 100644
--- a/Carpet/Carpet/src/Restrict.cc
+++ b/Carpet/Carpet/src/Restrict.cc
@@ -8,7 +8,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.5 2001/12/09 16:41:53 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.6 2002/01/09 17:45:39 schnetter Exp $";
 
 
 
@@ -24,36 +24,49 @@ namespace Carpet {
     
     Checkpoint ("%*sRestrict", 2*reflevel, "");
     
-    if (reflevel == hh->reflevels()-1) return;
-    
+    // Loop over variables with storage
     for (int group=0; group<CCTK_NumGroups(); ++group) {
-      
-      // Restrict only groups with storage
       if (CCTK_QueryGroupStorageI(cgh, group)) {
-	
-	const int tl = 0;
-	
 	for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
-	  for (int c=0; c<hh->components(reflevel); ++c) {
-	    arrdata[group].data[var]->ref_restrict (tl, reflevel, c, mglevel);
-	  }
-	  for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
-	    arrdata[group].data[var]->sync (tl, reflevel, c, mglevel);
-	  }
 	  
-	  for (int c=0; c<arrdata[group].hh->components(reflevel+1); ++c) {
-	    arrdata[group].data[var]->ref_bnd_prolongate
-	      (tl, reflevel+1, c, mglevel);
-	  }
-	  // TODO: is this necessary?
-	  for (int c=0; c<arrdata[group].hh->components(reflevel+1); ++c) {
-	    arrdata[group].data[var]->sync (tl, reflevel+1, c, mglevel);
-	  }
-	}
-	
-      }	// if group has storage
-      
+	  const int tl = 0;
+	  
+	  if (mglevel > 0) {
+	    
+	    for (int c=0; c<hh->components(reflevel); ++c) {
+	      arrdata[group].data[var]->mg_restrict (tl, reflevel, c, mglevel);
+	    }
+	    for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
+	      arrdata[group].data[var]->sync (tl, reflevel, c, mglevel);
+	    }
+	    
+	  } // if not finest multigrid level
+	  
+	  if (reflevel < hh->reflevels()-1) {
+	    
+	    for (int c=0; c<hh->components(reflevel); ++c) {
+	      arrdata[group].data[var]->ref_restrict
+		(tl, reflevel, c, mglevel);
+	    }
+	    for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
+	      arrdata[group].data[var]->sync (tl, reflevel, c, mglevel);
+	    }
+	    
+	    for (int c=0; c<arrdata[group].hh->components(reflevel+1); ++c) {
+	      arrdata[group].data[var]->ref_bnd_prolongate
+		(tl, reflevel+1, c, mglevel);
+	    }
+	    // TODO: is this necessary?
+	    for (int c=0; c<arrdata[group].hh->components(reflevel+1); ++c) {
+	      arrdata[group].data[var]->sync (tl, reflevel+1, c, mglevel);
+	    }
+	    
+	  } // if not finest refinement level
+	  
+	} // loop over variables
+      } // if group has storage
     } // loop over groups
+    
   }
   
 } // namespace Carpet
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 1cb27c766..faa8a4561 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -12,7 +12,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.18 2002/01/09 13:56:25 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.19 2002/01/09 17:45:40 schnetter Exp $";
 
 
 
@@ -46,6 +46,7 @@ namespace Carpet {
     // Multigrid information
     mglevels = multigrid_levels;
     mgfact = multigrid_factor;
+    maxmglevelfact = ipow(mgfact, mglevels-1);
     
     // Ghost zones
     vect<int,dim> lghosts, ughosts;
@@ -272,8 +273,8 @@ namespace Carpet {
     iteration.resize(maxreflevels, 0);
     
     // Set current position (this time for real)
-    set_reflevel  (cgh, 0);
     set_mglevel   (cgh, 0);
+    set_reflevel  (cgh, 0);
     set_component (cgh, -1);
     
     // Enable storage for all groups if desired
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index 5f097ae33..c16fd743a 100644
--- a/Carpet/Carpet/src/carpet_public.hh
+++ b/Carpet/Carpet/src/carpet_public.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.13 2002/01/09 13:56:25 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.14 2002/01/09 17:45:40 schnetter Exp $
 
 // It is assumed that the number of components of all arrays is equal
 // to the number of components of the grid functions, and that their
@@ -32,29 +32,35 @@ namespace Carpet {
   // Handle from CCTK_RegisterGHExtension
   extern int GHExtension;
   
+  // Multigrid levels
+  extern int mglevels;
+  
+  // Multigrid factor
+  extern int mgfact;
+  
   // Maximum number of refinement levels
   extern int maxreflevels;
   
+  // Multigrid factor on coarsest grid
+  extern int maxmglevelfact;
+  
   // Refinement factor
   extern int reffact;
   
   // Refinement factor on finest grid
   extern int maxreflevelfact;
   
-  // Multigrid levels
-  extern int mglevels;
-  
-  // Multigrid factor
-  extern int mgfact;
-  
   // Current iteration per refinement level
   extern vector<int> iteration;
   
   // Current position on the grid hierarchy
-  extern int reflevel;
   extern int mglevel;
+  extern int reflevel;
   extern int component;
   
+  // Current multigrid factor
+  extern int mglevelfact;
+  
   // Current refinement factor
   extern int reflevelfact;
   
@@ -148,6 +154,24 @@ namespace Carpet {
   
   
   
+  // Multigrid level iterator
+  
+#define BEGIN_MGLEVEL_LOOP(cgh)			\
+  do {						\
+    int _mgl;					\
+    assert (mglevel==0);			\
+    for (int _ml=mglevels-1; _ml>=0; --_ml) {	\
+      set_mglevel ((cGH*)(cgh), _ml);		\
+      {
+#define END_MGLEVEL_LOOP(cgh)			\
+      }						\
+    }						\
+    assert (mglevel==0);			\
+    _mgl = 0;					\
+  } while (0)
+  
+  
+  
   // Refinement level iterator
   
 #define BEGIN_REFLEVEL_LOOP(cgh)			\
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index 6209d1256..15a8c9f27 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -12,7 +12,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.15 2002/01/09 13:56:26 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.16 2002/01/09 17:45:40 schnetter Exp $";
 
 
 
@@ -176,47 +176,49 @@ namespace Carpet {
   
   
   
+  void set_mglevel (cGH* cgh, const int ml)
+  {
+    // Check
+    assert (ml>=0 && ml<mglevels);
+    assert (reflevel == 0);
+    assert (component == -1);
+    
+    mglevel = ml;
+    mglevelfact = ipow(mgfact, mglevel);
+    cgh->cctk_convlevel = mglevel;
+  }
+  
+  
+  
   void set_reflevel (cGH* cgh, const int rl)
   {
     // Check
-    assert (rl>=0 && rl<hh->reflevels());
+    assert (rl>=0 && rl<maxreflevels && rl<hh->reflevels());
     assert (component == -1);
     
     // Change
     reflevel = rl;
     const bbox<int,dim>& base = hh->baseextent;
-    reflevelfact = ipow(hh->reffact, reflevel);
-    cgh->cctk_delta_time = base_delta_time / reflevelfact;
-    for (int d=0; d<dim; ++d) {
-      cgh->cctk_gsh[d]
-	= ((base.shape() / base.stride()
-	    + dd->lghosts + dd->ughosts)[d] - 1) * reflevelfact + 1;
-      cgh->cctk_levfac[d] = reflevelfact;
-    }
+    reflevelfact = ipow(reffact, reflevel);
+    cgh->cctk_delta_time = base_delta_time / reflevelfact * mglevelfact;
+    assert (all(base.shape() % base.stride() == 0));
+    assert (all((base.shape() / base.stride()) % mglevelfact == 0));
+    vect<int,dim>::ref(cgh->cctk_gsh)
+      = (((base.shape() / base.stride() - 1) / mglevelfact
+	  + dd->lghosts + dd->ughosts)
+	 * reflevelfact + 1);
+    vect<int,dim>::ref(cgh->cctk_levfac) = reflevelfact;
     for (int group=0; group<CCTK_NumGroups(); ++group) {
       const bbox<int,dim>& base = arrdata[group].hh->baseextent;
-      for (int d=0; d<dim; ++d) {
-	((int*)arrdata[group].info.gsh)[d]
-	  = (((base.shape() / base.stride()
-	       + arrdata[group].dd->lghosts + arrdata[group].dd->ughosts)[d]
-	      - 1)
-	     * reflevelfact + 1);
-      }
+      vect<int,dim>::ref((int*)arrdata[group].info.gsh)
+	= (((base.shape() / base.stride() - 1) / mglevelfact
+	    + arrdata[group].dd->lghosts + arrdata[group].dd->ughosts)
+	   * reflevelfact + 1);
     }
   }
   
   
   
-  void set_mglevel (cGH* cgh, const int ml)
-  {
-    assert (ml==0);
-    assert (component==-1);
-    mglevel = ml;
-    cgh->cctk_convlevel = mglevel;
-  }
-  
-  
-  
   void set_component (cGH* cgh, const int c)
   {
     assert (reflevel>=0 && reflevel<hh->reflevels());
@@ -299,18 +301,28 @@ namespace Carpet {
 	assert (reflevel < (int)dd->boxes.size());
 	assert (component < (int)dd->boxes[reflevel].size());
 	assert (mglevel < (int)dd->boxes[reflevel][component].size());
+	const bbox<int,dim>& bext = hh->baseextent;
+	const bbox<int,dim>& iext = hh->extents[reflevel][component][mglevel];
 	const bbox<int,dim>& ext
 	  = dd->boxes[reflevel][component][mglevel].exterior;
 	for (int d=0; d<dim; ++d) {
 	  cgh->cctk_lsh[d] = (ext.shape() / ext.stride())[d];
 	  cgh->cctk_lbnd[d] = (ext.lower() / ext.stride())[d];
 	  cgh->cctk_ubnd[d] = (ext.upper() / ext.stride())[d];
-	  assert (cgh->cctk_lsh[d]>=0 && cgh->cctk_lsh[d]<=cgh->cctk_gsh[d]);
-	  assert (cgh->cctk_lbnd[d]>=0 && cgh->cctk_ubnd[d]<cgh->cctk_gsh[d]);
+ 	  assert (cgh->cctk_lsh[d]>=0 && cgh->cctk_lsh[d]<=cgh->cctk_gsh[d]);
+// 	  assert (cgh->cctk_lbnd[d]>=0 && cgh->cctk_ubnd[d]<cgh->cctk_gsh[d]);
+	  assert (cgh->cctk_ubnd[d]-cgh->cctk_lbnd[d]+1 == cgh->cctk_lsh[d]);
 	  assert (cgh->cctk_lbnd[d]<=cgh->cctk_ubnd[d]+1);
-	  // Do allow outer boundaries on the finer grids
-	  cgh->cctk_bbox[2*d  ] = cgh->cctk_lbnd[d] == 0;
-	  cgh->cctk_bbox[2*d+1] = cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1;
+	  // Do not allow outer boundaries on the finer grids
+	  if (reflevel==0) {
+	    assert (iext.lower()[d] >= bext.lower()[d]);
+	    assert (iext.upper()[d] <= bext.upper()[d]);
+	    cgh->cctk_bbox[2*d  ] = iext.lower()[d] == bext.lower()[d];
+	    cgh->cctk_bbox[2*d+1] = iext.upper()[d] == bext.upper()[d];
+	  } else {
+	    cgh->cctk_bbox[2*d  ] = 0;
+	    cgh->cctk_bbox[2*d+1] = 0;
+	  }
 	  for (int stg=0; stg<CCTK_NSTAGGER; ++stg) {
 	    // TODO: support staggering
 	    cgh->cctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d];
@@ -321,6 +333,9 @@ namespace Carpet {
 	assert (reflevel < (int)arrdata[group].dd->boxes.size());
 	assert (component < (int)arrdata[group].dd->boxes[reflevel].size());
 	assert (mglevel < (int)arrdata[group].dd->boxes[reflevel][component].size());
+	const bbox<int,dim>& bext = arrdata[group].hh->baseextent;
+	const bbox<int,dim>& iext
+	  = arrdata[group].hh->extents[reflevel][component][mglevel];
 	const bbox<int,dim>& ext
 	  = arrdata[group].dd->boxes[reflevel][component][mglevel].exterior;
 	for (int d=0; d<dim; ++d) {
@@ -332,14 +347,23 @@ namespace Carpet {
 	    = (ext.upper() / ext.stride())[d];
 	  assert (arrdata[group].info.lsh[d]>=0
 		  && arrdata[group].info.lsh[d]<=arrdata[group].info.gsh[d]);
-	  assert (arrdata[group].info.lbnd[d]>=0
-		  && arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]);
+// 	  assert (arrdata[group].info.lbnd[d]>=0
+// 		  && arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]);
+	  assert (arrdata[group].info.ubnd[d]-arrdata[group].info.lbnd[d]+1
+		  == arrdata[group].info.lsh[d]);
 	  assert (arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]+1);
-	  // Do allow outer boundaries on the finer grids
-	  ((int*)arrdata[group].info.bbox)[2*d  ]
-	    = arrdata[group].info.lbnd[d] == 0;
-	  ((int*)arrdata[group].info.bbox)[2*d+1]
-	    = arrdata[group].info.ubnd[d] == arrdata[group].info.gsh[d]-1;
+	  // Do not allow outer boundaries on the finer grids
+	  if (reflevel==0) {
+	    assert (iext.lower()[d] >= bext.lower()[d]);
+	    assert (iext.upper()[d] <= bext.upper()[d]);
+	    ((int*)arrdata[group].info.bbox)[2*d  ]
+	      = iext.lower()[d] == bext.lower()[d];
+	    ((int*)arrdata[group].info.bbox)[2*d+1]
+	      = iext.upper()[d] == bext.upper()[d];
+	  } else {
+	    ((int*)arrdata[group].info.bbox)[2*d  ] = 0;
+	    ((int*)arrdata[group].info.bbox)[2*d+1] = 0;
+	  }
 	}
       }
       
diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc
index a47073da2..a8a2a7c75 100644
--- a/Carpet/Carpet/src/variables.cc
+++ b/Carpet/Carpet/src/variables.cc
@@ -5,7 +5,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.4 2002/01/09 13:56:26 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.5 2002/01/09 17:45:40 schnetter Exp $";
 
 
 
@@ -18,6 +18,15 @@ namespace Carpet {
   // Handle from CCTK_RegisterGHExtension
   int GHExtension;
   
+  // Multigrid levels
+  int mglevels;
+  
+  // Multigrid factor
+  int mgfact;
+  
+  // Multigrid factor on coarsest grid
+  int maxmglevelfact;
+  
   // Maximum number of refinement levels
   int maxreflevels;
   
@@ -27,12 +36,6 @@ namespace Carpet {
   // Refinement factor on finest grid
   int maxreflevelfact;
   
-  // Multigrid levels
-  int mglevels;
-  
-  // Multigrid factor
-  int mgfact;
-  
   // Current iteration per refinement level
   vector<int> iteration;
   
@@ -41,6 +44,9 @@ namespace Carpet {
   int reflevel;
   int component;
   
+  // multigrid factor of current level: ipow(multigrid_factor, mglevel)
+  int mglevelfact;
+  
   // refinement factor of current level: ipow(refinement_factor, reflevel)
   int reflevelfact;
   
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
index cff34ddfb..62509615e 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.5 2001/12/09 16:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.6 2002/01/09 17:45:40 schnetter Exp $
       
 #include "cctk.h"
       
@@ -125,7 +125,7 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         if (fac.ne.0) then
-                           res = res + fac * src(i0+ii ,j0+jj ,k0+kk)
+                           res = res + fac * src(i0+ii, j0+jj, k0+kk)
                         end if
                         
                      end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
index 7b6e43da4..6cb18abbb 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.5 2002/01/09 17:45:41 schnetter Exp $
       
 #include "cctk.h"
       
@@ -142,8 +142,8 @@ c     Loop over fine region
                         
                         if (fac.ne.0) then
                            res = res
-     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s1fac * src1(i0+ii, j0+jj, k0+kk)
+     $                          + fac * s2fac * src2(i0+ii, j0+jj, k0+kk)
                         end if
                         
                      end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
index d6555ea0f..00a286952 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.6 2002/01/08 21:29:02 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.7 2002/01/09 17:45:41 schnetter Exp $
       
 #include "cctk.h"
       
@@ -157,8 +157,8 @@ c     Loop over fine region
                         
                         if (fac.ne.0) then
                            res = res
-     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1)
+     $                          + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1)
                         end if
                         
                      end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
index 0ec7eb6ce..f9c451070 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.2 2001/12/09 16:43:11 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.3 2002/01/09 17:45:41 schnetter Exp $
       
 #include "cctk.h"
       
@@ -37,9 +37,11 @@ c     bbox(:,3) is stride
       
       CCTK_REAL8 s1fac, s2fac, s3fac
       
+      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
       integer fi, fj, fk
+      integer ifac(2), jfac(2), kfac(2)
       integer ii, jj, kk
       integer fac
       CCTK_REAL8 res
@@ -113,55 +115,46 @@ c     Quadratic (second order) interpolation
       
       
 c     Loop over fine region
+      dstdiv = one / (dstifac * dstjfac * dstkfac)
+      
       do k = 0, regkext-1
+         k0 = (srckoff + k) / dstkfac
+         fk = mod(srckoff + k, dstkfac)
+         kfac(1) = (fk-dstkfac) * (-1)
+         kfac(2) = (fk        ) * 1
+         
          do j = 0, regjext-1
+            j0 = (srcjoff + j) / dstjfac
+            fj = mod(srcjoff + j, dstjfac)
+            jfac(1) = (fj-dstjfac) * (-1)
+            jfac(2) = (fj        ) * 1
+            
             do i = 0, regiext-1
-               
-               i0 = (srcioff + i) / dstifac + 1
-               j0 = (srcjoff + j) / dstjfac + 1
-               k0 = (srckoff + k) / dstkfac + 1
-               
+               i0 = (srcioff + i) / dstifac
                fi = mod(srcioff + i, dstifac)
-               fj = mod(srcjoff + j, dstjfac)
-               fk = mod(srckoff + k, dstkfac)
+               ifac(1) = (fi-dstifac) * (-1)
+               ifac(2) = (fi        ) * 1
                
                res = 0
                
-               do kk=0,1
-                  do jj=0,1
-                     do ii=0,1
-                        
-                        fac = 1
+               do kk=1,2
+                  do jj=1,2
+                     do ii=1,2
                         
-                        if (ii.eq.0) then
-                           fac = fac * (dstifac - fi)
-                        else
-                           fac = fac *            fi
-                        end if
-                        if (jj.eq.0) then
-                           fac = fac * (dstjfac - fj)
-                        else
-                           fac = fac *            fj
-                        end if
-                        if (kk.eq.0) then
-                           fac = fac * (dstkfac - fk)
-                        else
-                           fac = fac *            fk
-                        end if
+                        fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         if (fac.ne.0) then
                            res = res
-     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s3fac * src3(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s1fac * src1(i0+ii, j0+jj, k0+kk)
+     $                          + fac * s2fac * src2(i0+ii, j0+jj, k0+kk)
+     $                          + fac * s3fac * src3(i0+ii, j0+jj, k0+kk)
                         end if
                         
                      end do
                   end do
                end do
                
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
-     $              = res / (dstifac * dstjfac * dstkfac)
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
index 6317ccd00..2ce036fea 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.6 2002/01/08 21:29:03 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.7 2002/01/09 17:45:41 schnetter Exp $
       
 #include "cctk.h"
       
@@ -160,9 +160,9 @@ c     Loop over fine region
                         
                         if (fac.ne.0) then
                            res = res
-     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
-     $                          + fac * s3fac * src3(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1)
+     $                          + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1)
+     $                          + fac * s3fac * src3(i0+ii-1, j0+jj-1, k0+kk-1)
                         end if
                         
                      end do
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index 0eabc9a6e..f53038741 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.7 2002/01/08 12:03:56 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.8 2002/01/09 17:45:42 schnetter Exp $
 
  ***************************************************************************/
 
@@ -84,6 +84,10 @@ public:
     for (int d=0; d<D; ++d) elt[d]=(T)a[d];
   }
   
+  static vect& ref (T* const x) {
+    return *(vect*)x;
+  }
+  
   static vect dir (const int d) {
     vect r=(T)0;
     r[d]=1;
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 1f3b96b16..da3ea6edc 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -19,7 +19,7 @@
 #include "carpet.hh"
 #include "regrid.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.3 2002/01/09 13:56:28 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.4 2002/01/09 17:45:42 schnetter Exp $";
 
 
 
@@ -43,6 +43,9 @@ namespace CarpetRegrid {
     // Return if no regridding is desired
     if (regrid_every == -1) return 0;
     
+    // Return if this is not the main hierarchy
+    if (mglevel != 0) return 0;
+    
     // Return if this is the finest possible level
     if (reflevel == maxreflevels-1) return 0;