Pulled regridding out of the driver proper.

Pulled regridding out of the driver proper.
Without the new regridding thorn, Carpet will now default to a single
level, i.e. the base level only.  Use the thorn CarpetRegrid and set
the parameter CarpetRegrid::refinement_levels to get refinement.

darcs-hash:20011214153904-07bb3-b765f596812c0321ab0f898f93a9fbd192c49c09.gz

11 files changed, 264 insertions, 286 deletions

diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc
index 780f1b6b3..539c48bec 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.5 2001/11/02 17:51:15 schnetter Exp $";
+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 $";
 
 
 
@@ -52,7 +52,7 @@ namespace Carpet {
   {
     DECLARE_CCTK_PARAMETERS;
     
-    Checkpoint ("starting Evolve...");
+    Waypoint ("starting Evolve...");
     
     const int convlev = 0;
     cGH* cgh = fc->GH[convlev];
@@ -64,7 +64,7 @@ namespace Carpet {
       ++cgh->cctk_iteration;
       cgh->cctk_time += base_delta_time / maxreflevelfact;
       
-      Checkpoint ("Evolving iteration %d...", cgh->cctk_iteration);
+      Waypoint ("Evolving iteration %d...", cgh->cctk_iteration);
       
       BEGIN_REFLEVEL_LOOP(cgh) {
 	if ((cgh->cctk_iteration-1) % (maxreflevelfact/reflevelfact) == 0) {
@@ -94,16 +94,19 @@ namespace Carpet {
 	  Poison (cgh, currenttimebutnotifonly);
 	  
 	  // Evolve
-	  Checkpoint ("%*sScheduling PRESTEP", 2*reflevel, "");
+	  Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
 	  CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
-	  Checkpoint ("%*sScheduling EVOL", 2*reflevel, "");
+	  Waypoint ("%*sScheduling EVOL", 2*reflevel, "");
 	  CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
-	  Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
+	  Waypoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
 	  CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
 	  
 	  // Checking
 	  PoisonCheck (cgh, currenttimebutnotifonly);
 	  
+	  // Recompose grid hierarchy
+	  Recompose (cgh);
+	  
 	}
       } END_REFLEVEL_LOOP(cgh);
       
@@ -116,16 +119,16 @@ namespace Carpet {
 	  // Checking
 	  CalculateChecksums (cgh, currenttime);
 	  
-	  // Checkpoint
-	  Checkpoint ("%*sScheduling CHECKPOINT", 2*reflevel, "");
+	  // Waypoint
+	  Waypoint ("%*sScheduling CHECKPOINT", 2*reflevel, "");
 	  CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction);
 	  
 	  // Analysis
-	  Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
+	  Waypoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
 	  CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
 	  
 	  // Output
-	  Checkpoint ("%*sOutputGH", 2*reflevel, "");
+	  Waypoint ("%*sOutputGH", 2*reflevel, "");
 	  CCTK_OutputGH (cgh);
 	  
 	  // Checking
@@ -136,7 +139,7 @@ namespace Carpet {
       
     } // main loop
     
-    Checkpoint ("done with Evolve.");
+    Waypoint ("done with Evolve.");
     
     return 0;
   }
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 45c6af58d..ffee2a04e 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.1 2001/07/04 12:29:47 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.2 2001/12/14 16:39:07 schnetter Exp $";
 
 
 
@@ -29,7 +29,7 @@ namespace Carpet {
     CCTKi_AddGH (fc, convlev, cgh);
     
     // Delay checkpoint until MPI has been initialised
-    Checkpoint ("starting Initialise...");
+    Waypoint ("starting Initialise...");
     
     // Initialise stuff
     cgh->cctk_iteration = 0;
@@ -42,17 +42,17 @@ namespace Carpet {
     CCTKi_InitGHExtensions (cgh);
     
     // Check parameters
-    Checkpoint ("PARAMCHECK");
+    Waypoint ("PARAMCHECK");
     CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction);
     CCTKi_FinaliseParamWarn();
     
     BEGIN_REFLEVEL_LOOP(cgh) {
       
       // Checking
-      Poison (cgh, allbutlasttime);
+      Poison (cgh, alltimes);
       
       // Set up the grid
-      Checkpoint ("%*sScheduling BASEGRID", 2*reflevel, "");
+      Waypoint ("%*sScheduling BASEGRID", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction);
       if (reflevel==0) {
 	base_delta_time = cgh->cctk_delta_time;
@@ -64,23 +64,26 @@ namespace Carpet {
       }
       
       // Set up the initial data
-      Checkpoint ("%*sScheduling INITIAL", 2*reflevel, "");
+      Waypoint ("%*sScheduling INITIAL", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction);
-      Checkpoint ("%*sScheduling POSTINITIAL", 2*reflevel, "");
+      Waypoint ("%*sScheduling POSTINITIAL", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction);
       
       // Recover
-      Checkpoint ("%*sScheduling RECOVER_VARIABLES", 2*reflevel, "");
+      Waypoint ("%*sScheduling RECOVER_VARIABLES", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_RECOVER_VARIABLES", cgh, CallFunction);
-      Checkpoint ("%*sScheduling CPINITIAL", 2*reflevel, "");
+      Waypoint ("%*sScheduling CPINITIAL", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction);
       
       // Poststep
-      Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
+      Waypoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
       
       // Checking
-      PoisonCheck (cgh, allbutlasttime);
+      PoisonCheck (cgh, alltimes);
+      
+      // Recompose grid hierarchy
+      Recompose (cgh);
       
     } END_REFLEVEL_LOOP(cgh);
     
@@ -90,22 +93,22 @@ namespace Carpet {
       Restrict (cgh);
       
       // Checking
-      CalculateChecksums (cgh, allbutlasttime);
+      CalculateChecksums (cgh, allbutcurrenttime);
       
       // Analysis
-      Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
+      Waypoint ("%*sScheduling ANALYSIS", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
       
       // Output
-      Checkpoint ("%*sOutputGH", 2*reflevel, "");
+      Waypoint ("%*sOutputGH", 2*reflevel, "");
       CCTK_OutputGH (cgh);
       
       // Checking
-      CheckChecksums (cgh, allbutlasttime);
+      CheckChecksums (cgh, allbutcurrenttime);
       
     } END_REVERSE_REFLEVEL_LOOP(cgh);
     
-    Checkpoint ("done with Initialise.");
+    Waypoint ("done with Initialise.");
     
     return 0;
   }
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index d9f2525f5..1fded5cb8 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -2,18 +2,20 @@
 #include <math.h>
 #include <stdlib.h>
 
+#include <list>
 #include <vector>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 
 #include "Carpet/CarpetLib/src/bbox.hh"
+#include "Carpet/CarpetLib/src/bboxset.hh"
 #include "Carpet/CarpetLib/src/gh.hh"
 #include "Carpet/CarpetLib/src/vect.hh"
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.11 2001/12/07 18:24:16 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.12 2001/12/14 16:39:07 schnetter Exp $";
 
 
 
@@ -29,15 +31,51 @@ namespace Carpet {
   
   
   
+  static void SplitRegions_AlongZ       (const cGH* cgh, gh<dim>::cexts& bbss);
+  static void SplitRegions_AsSpecified  (const cGH* cgh, gh<dim>::cexts& bbss);
+  
+  static void MakeProcessors_RoundRobin (const cGH* cgh,
+					 const gh<dim>::rexts& bbss,
+					 gh<dim>::rprocs& pss);
+  
+  
+  
+  static bool RegionsOK (const gh<dim>::rexts& bbsss,
+			 const gh<dim>::rprocs& pss);
   static void Adapt (const cGH* cgh, int reflevels, gh<dim>* hh);
   static void Output (const cGH* cgh, const gh<dim>* hh, const char* descr);
   
   
   
+  bool RegionsOK (const gh<dim>::rexts& bbsss, const gh<dim>::rprocs& pss)
+  {
+    // At least one level
+    if (bbsss.size() == 0) return false;
+    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
+      // No empty levels
+      if (bbsss[rl].size() == 0) return false;
+      for (int c=0; c<(int)bbsss[rl].size(); ++c) {
+	// At least one multigrid level
+	if (bbsss[rl][c].size() == 0) return false;
+      }
+    }
+    
+    if (pss.size() != bbsss.size()) return false;
+    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
+      if (pss[rl].size() != bbsss[rl].size()) return false;
+    }
+    
+    return true;
+  }
+  
+  
+  
   void RegisterRecomposeRegions (const gh<dim>::rexts& bbsss,
 				 const gh<dim>::rprocs& pss)
   {
-    // save the region information for the next regridding
+    // Check the regions
+    assert (RegionsOK(bbsss,pss));
+    // Save the region information for the next regridding
     next_bbsss = bbsss;
     next_pss = pss;
     do_recompose = true;
@@ -50,11 +88,14 @@ namespace Carpet {
     DECLARE_CCTK_PARAMETERS;
     
     assert (component == -1);
-    Checkpoint ("%*sRecompose", 2*reflevel, "");
+    Waypoint ("%*sRecompose", 2*reflevel, "");
     
     // Check whether to recompose
     if (!do_recompose) return;
     
+    // Check the regions
+    assert (RegionsOK(next_bbsss, next_pss));
+    
     // Recompose
     hh->recompose (next_bbsss, next_pss);
     Output (cgh, hh, 0);
@@ -153,226 +194,172 @@ namespace Carpet {
   
   
   
-  // This is a helpful helper routine. The user can use it to define
-  // how the hierarchy should be refined. But the result of this
-  // routine is rather arbitrary.
-  void MakeRegions_RefineCentre (const cGH* cgh, const int reflevels,
-				 gh<dim>::rexts& bbsss)
+  static void Output (const cGH* cgh, const gh<dim>* hh, const char* descr)
   {
     DECLARE_CCTK_PARAMETERS;
     
-    const int mglevels = 1;	// arbitrary value
-    
-    vector<vector<bbox<int,dim> > > bbss(reflevels);
-    // note: what this routine calls "ub" is "ub+str" elsewhere
-    vect<int,dim> rstr = hh->baseextent.stride();
-    vect<int,dim> rlb  = hh->baseextent.lower();
-    vect<int,dim> rub  = hh->baseextent.upper() + rstr;
-    
-    for (int rl=0; rl<reflevels; ++rl) {
-      if (rl>0) {
-	// save old values
-	const vect<int,dim> oldrlb = rlb;
-	const vect<int,dim> oldrub = rub;
-	// calculate extent and centre
-	const vect<int,dim> rextent = rub - rlb;
-	const vect<int,dim> rcentre = rlb + (rextent / 2 / rstr) * rstr;
-	// calculate new extent
-	assert (all(rextent % hh->reffact == 0));
-	const vect<int,dim> newrextent = rextent / hh->reffact;
-	// refined boxes have smaller stride
-	assert (all(rstr%hh->reffact == 0));
-	rstr /= hh->reffact;
-	// refine (arbitrarily) around the center only
-	rlb = rcentre - (newrextent/2 / rstr) * rstr;
-// 	// refine (arbitrarily) around the lower boundary only
-// 	rlb = rlb;
-	rub = rlb + newrextent;
-	// require rub<oldrub because we really want rub-rstr<=oldrub-oldstr
-	assert (all(rlb >= oldrlb && rub < oldrub));
+    if (verbose && CCTK_MyProc(cgh)==0) {
+      cout << endl;
+      cout << "New bounding boxes";
+      if (descr) {
+	if (strlen(descr)) {
+	  cout << " for group " << descr;
+	} else {
+	  cout << " for scalars";
+	}
       }
-      vector<bbox<int,dim> > bbs(1);
-      bbs[0] = bbox<int,dim>(rlb, rub-rstr, rstr);
-      bbss[rl] = bbs;
+      cout  << ":" << endl;
+      for (int rl=0; rl<hh->reflevels(); ++rl) {
+	for (int c=0; c<hh->components(rl); ++c) {
+	  for (int ml=0; ml<hh->mglevels(rl,c); ++ml) {
+	    cout << "   rl " << rl << "   c " << c << "   ml " << ml
+		 << "   bbox " << hh->extents[rl][c][ml] << endl;
+	  }
+	}
+      }
+      cout << endl;
+      cout << "New processor distribution";
+      if (descr) {
+	if (strlen(descr)) {
+	  cout << " for group " << descr;
+	} else {
+	  cout << " for scalars";
+	}
+      }
+      cout  << ":" << endl;
+      for (int rl=0; rl<hh->reflevels(); ++rl) {
+	for (int c=0; c<hh->components(rl); ++c) {
+	  cout << "   rl " << rl << "   c " << c
+	       << "   processor " << hh->processors[rl][c] << endl;
+	}
+      }
+      cout << endl;
     }
-    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
   }
   
   
   
-  void MakeRegions_AsSpecified (const cGH* cgh, const int reflevels,
-				gh<dim>::rexts& bbsss)
+  void SplitRegions (const cGH* cgh, gh<dim>::cexts& bbss)
   {
     DECLARE_CCTK_PARAMETERS;
     
-    const int mglevels = 1;	// arbitrary value
-    
-    const vect<int,dim> rstr = hh->baseextent.stride();
-    const vect<int,dim> rlb  = hh->baseextent.lower();
-    const vect<int,dim> rub  = hh->baseextent.upper();
-    
-    if (reflevels>4) {
-      CCTK_WARN (0, "Cannot currently specify refinement regions for more than 4 refinement levels");
-    }
-    
-    assert (reflevels<4);
-    vector<vect<int,dim> > lower(4), upper(4);
-    lower[0] = rlb;
-    upper[0] = rub;
-    lower[1] = vect<int,dim> (l1xmin, l1ymin, l1zmin);
-    upper[1] = vect<int,dim> (l1xmax, l1ymax, l1zmax);
-    lower[2] = vect<int,dim> (l2xmin, l2ymin, l2zmin);
-    upper[2] = vect<int,dim> (l2xmax, l2ymax, l2zmax);
-    lower[3] = vect<int,dim> (l3xmin, l3ymin, l3zmin);
-    upper[3] = vect<int,dim> (l3xmax, l3ymax, l3zmax);
-    
-    vector<vector<bbox<int,dim> > > bbss(reflevels);
-    
-    for (int rl=0; rl<reflevels; ++rl) {
-      const int levfac = floor(pow((double)hh->reffact, rl) + 0.5);
-      assert (all (rstr % levfac == 0));
-      const vect<int,dim> str (rstr / levfac);
-      const vect<int,dim> lb  (lower[rl]);
-      const vect<int,dim> ub  (upper[rl]);
-      if (! all(lb>=rlb && ub<=rub)) {
-	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "The refinement region boundaries for refinement level #%d are not within the main grid", rl);
-      }
-      if (! all(lb<=ub)) {
-	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "The refinement region boundaries for refinement level #%d have the upper boundary less than the lower boundary", rl);
-      }
-      if (! all(lb%str==0 && ub%str==0)) {
-	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
-      }
-      assert (all(lb>=rlb && ub<=rub));
-      assert (all(lb<=ub));
-      assert (all(lb%str==0 && ub%str==0));
-      vector<bbox<int,dim> > bbs(1);
-      bbs[0] = bbox<int,dim>(lb, ub, str);
-      bbss[rl] = bbs;
-    }
-    
-    if (reflevels>0) {
-      assert (bbss[0].size() > 0);
-      assert (bbss[0][0] == hh->baseextent);
+    if (CCTK_EQUALS (processor_topology, "automatic")) {
+      SplitRegions_AlongZ (cgh, bbss);
+    } else if (CCTK_EQUALS (processor_topology, "manual")) {
+      SplitRegions_AsSpecified (cgh, bbss);
+    } else {
+      abort();
     }
-    
-    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
   }
   
   
   
-  // This is a helpful helper routine. The user can use it to define
-  // how the hierarchy should be refined. But the result of this
-  // routine is rather arbitrary.
-  void SplitRegions_AlongZ (const cGH* cgh, gh<dim>::rexts& bbsss)
+  void SplitRegions_AlongZ (const cGH* cgh, gh<dim>::cexts& bbss)
   {
-    DECLARE_CCTK_PARAMETERS;
+    // Something to do?
+    if (bbss.size() == 0) return;
     
     const int nprocs    = CCTK_nProcs(cgh);
     const int mglevels  = 1;	// arbitrary value
     
-    vector<vector<bbox<int,dim> > > bbss(bbsss.size());
+    assert (bbss.size() == 1);
+    assert (bbss[0].size() == 1);
     
-    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
-      assert (bbsss[rl].size() == 1);
-      assert (bbsss[rl][0].size() == 1);
-      
-      const vect<int,dim> rstr = bbsss[rl][0][0].stride();
-      const vect<int,dim> rlb  = bbsss[rl][0][0].lower();
-      const vect<int,dim> rub  = bbsss[rl][0][0].upper() + rstr;
-      
-      vector<bbox<int,dim> > bbs(nprocs);
-      for (int c=0; c<nprocs; ++c) {
-	vect<int,dim> cstr = rstr;
-	vect<int,dim> clb = rlb;
-	vect<int,dim> cub = rub;
-	const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1];
-	const int locnpz = (glonpz + nprocs - 1) / nprocs;
-	const int zstep = locnpz * cstr[dim-1];
-	clb[dim-1] = rlb[dim-1] + zstep *  c;
-	cub[dim-1] = rlb[dim-1] + zstep * (c+1);
- 	if (clb[dim-1] > rub[dim-1]) clb[dim-1] = rub[dim-1];
- 	if (cub[dim-1] > rub[dim-1]) cub[dim-1] = rub[dim-1];
-	assert (clb[dim-1] <= cub[dim-1]);
-	assert (cub[dim-1] <= rub[dim-1]);
-	bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
-      }
-      bbss[rl] = bbs;
+    const vect<int,dim> rstr = bbss[0][0].stride();
+    const vect<int,dim> rlb  = bbss[0][0].lower();
+    const vect<int,dim> rub  = bbss[0][0].upper() + rstr;
+    
+    vector<bbox<int,dim> > bbs(nprocs);
+    for (int c=0; c<nprocs; ++c) {
+      vect<int,dim> cstr = rstr;
+      vect<int,dim> clb = rlb;
+      vect<int,dim> cub = rub;
+      const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1];
+      const int locnpz = (glonpz + nprocs - 1) / nprocs;
+      const int zstep = locnpz * cstr[dim-1];
+      clb[dim-1] = rlb[dim-1] + zstep *  c;
+      cub[dim-1] = rlb[dim-1] + zstep * (c+1);
+      if (clb[dim-1] > rub[dim-1]) clb[dim-1] = rub[dim-1];
+      if (cub[dim-1] > rub[dim-1]) cub[dim-1] = rub[dim-1];
+      assert (clb[dim-1] <= cub[dim-1]);
+      assert (cub[dim-1] <= rub[dim-1]);
+      bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
     }
-    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
+    bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
   }
   
   
   
-  void SplitRegions_AsSpecified (const cGH* cgh, gh<dim>::rexts& bbsss)
+  void SplitRegions_AsSpecified (const cGH* cgh, gh<dim>::cexts& bbss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    // Something to do?
+    if (bbss.size() == 0) return;
+    
     const int nprocs    = CCTK_nProcs(cgh);
     const int mglevels  = 1;	// arbitrary value
     
-    vector<vector<bbox<int,dim> > > bbss(bbsss.size());
+    assert (bbss.size() == 1);
+    assert (bbss[0].size() == 1);
     
-    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
-      assert (bbsss[rl].size() == 1);
-      assert (bbsss[rl][0].size() == 1);
-      
-      const vect<int,dim> rstr = bbsss[rl][0][0].stride();
-      const vect<int,dim> rlb  = bbsss[rl][0][0].lower();
-      const vect<int,dim> rub  = bbsss[rl][0][0].upper() + rstr;
-      
-      const vect<int,dim> nprocs_dir
-	(processor_topology_3d_x, processor_topology_3d_y,
-	 processor_topology_3d_z);
-      assert (all (nprocs_dir > 0));
-      if (prod(nprocs_dir) != nprocs) {
-	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "The specified processor topology [%d,%d,%d] does not fit the number of processors, which is %d", nprocs_dir[0], nprocs_dir[1], nprocs_dir[2], nprocs);
-      }
-      assert (prod(nprocs_dir) == nprocs);
-      
-      vector<bbox<int,dim> > bbs(nprocs);
-      assert (dim==3);
-      for (int k=0; k<nprocs_dir[2]; ++k) {
-	for (int j=0; j<nprocs_dir[1]; ++j) {
-	  for (int i=0; i<nprocs_dir[0]; ++i) {
-	    const int c = i + nprocs_dir[0] * (j + nprocs_dir[1] * k);
-	    const vect<int,dim> ipos (i, j, k);
-	    vect<int,dim> cstr = rstr;
-	    vect<int,dim> clb = rlb;
-	    vect<int,dim> cub = rub;
-	    const vect<int,dim> glonp = (rub - rlb) / cstr;
-	    const vect<int,dim> locnp = (glonp + nprocs_dir - 1) / nprocs_dir;
-	    const vect<int,dim> step = locnp * cstr;
-	    clb = rlb + step *  ipos;
-	    cub = rlb + step * (ipos+1);
-	    clb = min (clb, rub);
-	    cub = min (cub, rub);
-	    assert (all (clb <= cub));
-	    assert (all (cub <= rub));
-	    bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
-	  }
+    const vect<int,dim> rstr = bbss[0][0].stride();
+    const vect<int,dim> rlb  = bbss[0][0].lower();
+    const vect<int,dim> rub  = bbss[0][0].upper() + rstr;
+    
+    const vect<int,dim> nprocs_dir
+      (processor_topology_3d_x, processor_topology_3d_y,
+       processor_topology_3d_z);
+    assert (all (nprocs_dir > 0));
+    if (prod(nprocs_dir) != nprocs) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "The specified processor topology [%d,%d,%d] does not fit the number of processors, which is %d", nprocs_dir[0], nprocs_dir[1], nprocs_dir[2], nprocs);
+    }
+    assert (prod(nprocs_dir) == nprocs);
+    
+    vector<bbox<int,dim> > bbs(nprocs);
+    assert (dim==3);
+    for (int k=0; k<nprocs_dir[2]; ++k) {
+      for (int j=0; j<nprocs_dir[1]; ++j) {
+	for (int i=0; i<nprocs_dir[0]; ++i) {
+	  const int c = i + nprocs_dir[0] * (j + nprocs_dir[1] * k);
+	  const vect<int,dim> ipos (i, j, k);
+	  vect<int,dim> cstr = rstr;
+	  vect<int,dim> clb = rlb;
+	  vect<int,dim> cub = rub;
+	  const vect<int,dim> glonp = (rub - rlb) / cstr;
+	  const vect<int,dim> locnp = (glonp + nprocs_dir - 1) / nprocs_dir;
+	  const vect<int,dim> step = locnp * cstr;
+	  clb = rlb + step *  ipos;
+	  cub = rlb + step * (ipos+1);
+	  clb = min (clb, rub);
+	  cub = min (cub, rub);
+	  assert (all (clb <= cub));
+	  assert (all (cub <= rub));
+	  bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
 	}
       }
-      bbss[rl] = bbs;
     }
-    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
+    bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
+  }
+  
+  
+  
+  void MakeProcessors (const cGH* cgh, const gh<dim>::rexts& bbsss,
+		       gh<dim>::rprocs& pss)
+  {
+    MakeProcessors_RoundRobin (cgh, bbsss, pss);
   }
   
   
   
   // This is a helpful helper routine. The user can use it to define
-  // how the hierarchy should be refined. But the result of this
+  // how the hierarchy should be refined.  But the result of this
   // routine is rather arbitrary.
   void MakeProcessors_RoundRobin (const cGH* cgh, const gh<dim>::rexts& bbsss,
 				  gh<dim>::rprocs& pss)
   {
-    DECLARE_CCTK_PARAMETERS;
-    
-    const int nprocs    = CCTK_nProcs(cgh);
+    const int nprocs = CCTK_nProcs(cgh);
     
     pss.resize(bbsss.size());
     for (int rl=0; rl<(int)bbsss.size(); ++rl) {
@@ -385,49 +372,4 @@ namespace Carpet {
     }
   }
   
-  
-  
-  static void Output (const cGH* cgh, const gh<dim>* hh, const char* descr)
-  {
-    DECLARE_CCTK_PARAMETERS;
-    
-    if (verbose && CCTK_MyProc(cgh)==0) {
-      cout << endl;
-      cout << "New bounding boxes";
-      if (descr) {
-	if (strlen(descr)) {
-	  cout << " for group " << descr;
-	} else {
-	  cout << " for scalars";
-	}
-      }
-      cout  << ":" << endl;
-      for (int rl=0; rl<hh->reflevels(); ++rl) {
-	for (int c=0; c<hh->components(rl); ++c) {
-	  for (int ml=0; ml<hh->mglevels(rl,c); ++ml) {
-	    cout << "   rl " << rl << "   c " << c << "   ml " << ml
-		 << "   bbox " << hh->extents[rl][c][ml] << endl;
-	  }
-	}
-      }
-      cout << endl;
-      cout << "New processor distribution";
-      if (descr) {
-	if (strlen(descr)) {
-	  cout << " for group " << descr;
-	} else {
-	  cout << " for scalars";
-	}
-      }
-      cout  << ":" << endl;
-      for (int rl=0; rl<hh->reflevels(); ++rl) {
-	for (int c=0; c<hh->components(rl); ++c) {
-	  cout << "   rl " << rl << "   c " << c
-	       << "   processor " << hh->processors[rl][c] << endl;
-	}
-      }
-      cout << endl;
-    }
-  }
-  
 } // namespace Carpet
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 2df49033c..1efb5e3b3 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -11,7 +11,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.15 2001/12/09 16:41:53 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.16 2001/12/14 16:39:08 schnetter Exp $";
 
 
 
@@ -35,7 +35,7 @@ namespace Carpet {
     CCTK_VInfo (CCTK_THORNSTRING,
 		"Carpet is running on %d processors", CCTK_nProcs(cgh));
     
-    Checkpoint ("starting SetupGH...");
+    Waypoint ("starting SetupGH...");
     
     // Refinement information
     maxreflevels = max_refinement_levels;
@@ -240,22 +240,20 @@ namespace Carpet {
     component = -1;
     
     // Invent a refinement structure
+    vector<bbox<int,dim> > bbs(1);
+    bbs[0] = hh->baseextent;
+    
+    vector<vector<bbox<int,dim> > > bbss(1);
+    bbss[0] = bbs;
+    
+    SplitRegions (cgh, bbss);
+    
     gh<dim>::rexts bbsss;
+    bbsss = hh->make_multigrid_boxes(bbss, 1);
+    
     gh<dim>::rprocs pss;
-    assert (refinement_levels <= max_refinement_levels);
-    if (CCTK_EQUALS (refined_regions, "centre")) {
-      MakeRegions_RefineCentre (cgh, refinement_levels, bbsss);
-    } else if (CCTK_EQUALS (refined_regions, "manual")) {
-      MakeRegions_AsSpecified (cgh, refinement_levels, bbsss);
-    }
-    if (CCTK_EQUALS (processor_topology, "automatic")) {
-      SplitRegions_AlongZ (cgh, bbsss);
-    } else if (CCTK_EQUALS (processor_topology, "manual")) {
-      SplitRegions_AsSpecified (cgh, bbsss);
-    } else {
-      abort();
-    }
-    MakeProcessors_RoundRobin (cgh, bbsss, pss);
+    MakeProcessors (cgh, bbsss, pss);
+    
     RegisterRecomposeRegions (bbsss, pss);
     
     // Recompose grid hierarchy
@@ -280,7 +278,7 @@ namespace Carpet {
       }
     }
     
-    Checkpoint ("done with SetupGH.");
+    Waypoint ("done with SetupGH.");
     
     // We register only once, ergo we get only one handle, ergo there
     // is only one grid hierarchy for us.  We store that statically,
diff --git a/Carpet/Carpet/src/Shutdown.cc b/Carpet/Carpet/src/Shutdown.cc
index 1e31d7919..b768b5e30 100644
--- a/Carpet/Carpet/src/Shutdown.cc
+++ b/Carpet/Carpet/src/Shutdown.cc
@@ -9,7 +9,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.3 2001/11/02 17:51:15 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.4 2001/12/14 16:39:08 schnetter Exp $";
 
 
 
@@ -23,20 +23,20 @@ namespace Carpet {
   {
     DECLARE_CCTK_PARAMETERS;
     
-    Checkpoint ("starting Shutdown...");
+    Waypoint ("starting Shutdown...");
     
     const int convlev = 0;
     cGH* cgh = fc->GH[convlev];
     
     // Terminate
     BEGIN_REFLEVEL_LOOP(cgh) {
-      Checkpoint ("%*sScheduling TERMINATE", 2*reflevel, "");
+      Waypoint ("%*sScheduling TERMINATE", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction);
     } END_REFLEVEL_LOOP(cgh);
     
     // Shutdown
     BEGIN_REFLEVEL_LOOP(cgh) {
-      Checkpoint ("%*sScheduling SHUTDOWN", 2*reflevel, "");
+      Waypoint ("%*sScheduling SHUTDOWN", 2*reflevel, "");
       CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction);
     } END_REFLEVEL_LOOP(cgh);
     
@@ -44,7 +44,7 @@ namespace Carpet {
     printf ("Done.\n");
     
     // earlier checkpoint before calling finalising MPI
-    Checkpoint ("done with Shutdown.");
+    Waypoint ("done with Shutdown.");
     
     dist::finalize();
     
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index b825b470c..18b49cc74 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -9,7 +9,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.6 2001/12/09 16:41:53 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.7 2001/12/14 16:39:08 schnetter Exp $";
 
 
 
@@ -55,6 +55,15 @@ namespace Carpet {
     const int my_prolongation_order_time
       = num_tl==1 ? 0 : prolongation_order_time;
     
+    if (max_refinement_levels > 1) {
+      if (num_tl <= my_prolongation_order_time) {
+	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+		    "There are not enough time levels for the desired temporal prolongation order in the grid function group \"%s\".  With Carpet::prolongation_order_time=%d, you need at least %d time levels.",
+		    CCTK_GroupName(group),
+		    prolongation_order_time, prolongation_order_time+1);
+      }
+    }
+    
     assert (arrdata[group].data.size()==0
 	    || arrdata[group].data[0] == 0);
     for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh
index 711091bc8..00ef19dfe 100644
--- a/Carpet/Carpet/src/carpet.hh
+++ b/Carpet/Carpet/src/carpet.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.13 2001/11/02 10:59:00 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.14 2001/12/14 16:39:09 schnetter Exp $
 
 #include "carpet_public.hh"
 
@@ -25,6 +25,7 @@ namespace Carpet {
   void CheckChecksums (cGH* cgh, checktimes where);
   
   // Debugging output
+  void Waypoint (const char* fmt, ...);
   void Checkpoint (const char* fmt, ...);
   
   // Error output
diff --git a/Carpet/Carpet/src/carpet_public.h b/Carpet/Carpet/src/carpet_public.h
index 38a494416..1fa8e4fe1 100644
--- a/Carpet/Carpet/src/carpet_public.h
+++ b/Carpet/Carpet/src/carpet_public.h
@@ -1,4 +1,7 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.h,v 1.4 2001/12/09 16:41:53 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.h,v 1.5 2001/12/14 16:39:09 schnetter Exp $ */
+
+#ifndef CARPET_PUBLIC_H
+#define CARPET_PUBLIC_H
 
 #include <mpi.h>
 
@@ -14,3 +17,5 @@ int CarpetStartup (void);
 /* Helper functions */
 MPI_Comm CarpetMPIComm (void);
 MPI_Datatype CarpetMPIDatatype (int vartype);
+
+#endif /* ! defined(CARPET_PUBLIC_H) */
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index fd0aa7a7a..bcb483926 100644
--- a/Carpet/Carpet/src/carpet_public.hh
+++ b/Carpet/Carpet/src/carpet_public.hh
@@ -1,16 +1,20 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.8 2001/12/09 16:41:53 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.9 2001/12/14 16:39:09 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
 // distribution onto the processors is the same, and that all
 // processors own the same number of components.
 
+#ifndef CARPET_PUBLIC_HH
+#define CARPET_PUBLIC_HH
+
 #include <vector>
 
 #include "cctk.h"
 #include "cctk_Schedule.h"
 
 #include "Carpet/CarpetLib/src/dh.hh"
+#include "Carpet/CarpetLib/src/gf.hh"
 #include "Carpet/CarpetLib/src/ggf.hh"
 #include "Carpet/CarpetLib/src/gh.hh"
 #include "Carpet/CarpetLib/src/th.hh"
@@ -121,14 +125,10 @@ namespace Carpet {
   void RegisterRecomposeRegions (const gh<dim>::rexts& bbsss,
 				 const gh<dim>::rprocs& pss);
   
-  void MakeRegions_RefineCentre  (const cGH* cgh, int reflevels,
-				  gh<dim>::rexts& bbsss);
-  void MakeRegions_AsSpecified   (const cGH* cgh, int reflevels,
-				  gh<dim>::rexts& bbsss);
-  void SplitRegions_AlongZ       (const cGH* cgh, gh<dim>::rexts& bbsss);
-  void SplitRegions_AsSpecified  (const cGH* cgh, gh<dim>::rexts& bbsss);
-  void MakeProcessors_RoundRobin (const cGH* cgh, const gh<dim>::rexts& bbsss,
-				  gh<dim>::rprocs& pss);
+  void SplitRegions (const cGH* cgh, gh<dim>::cexts& bbss);
+  
+  void MakeProcessors (const cGH* cgh, const gh<dim>::rexts& bbsss,
+		       gh<dim>::rprocs& pss);
   
   
   
@@ -200,3 +200,5 @@ namespace Carpet {
   } while (0)
   
 } // namespace Carpet
+
+#endif // ! defined(CARPET_PUBLIC_HH)
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index b2b6bc88e..fb36ba70c 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -11,7 +11,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.12 2001/12/09 16:41:53 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.13 2001/12/14 16:39:10 schnetter Exp $";
 
 
 
@@ -129,10 +129,26 @@ namespace Carpet {
   
   
   
+  void Waypoint (const char* fmt, ...)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    if (verbose || veryverbose) {
+      va_list args;
+      char msg[1000];
+      va_start (args, fmt);
+      vsnprintf (msg, sizeof(msg), fmt, args);
+      va_end (args);
+      CCTK_INFO (msg);
+    }
+    if (barriers) {
+      MPI_Barrier (dist::comm);
+    }
+  }
+  
   void Checkpoint (const char* fmt, ...)
   {
     DECLARE_CCTK_PARAMETERS;
-    if (verbose) {
+    if (veryverbose) {
       va_list args;
       char msg[1000];
       va_start (args, fmt);
diff --git a/Carpet/Carpet/src/make.code.defn b/Carpet/Carpet/src/make.code.defn
index 446f04cc6..8086ed933 100644
--- a/Carpet/Carpet/src/make.code.defn
+++ b/Carpet/Carpet/src/make.code.defn
@@ -1,9 +1,8 @@
 # Main make.code.defn file for thorn Carpet   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.3 2001/12/05 03:31:57 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.4 2001/12/14 16:39:10 schnetter Exp $
 
 # Source files in this directory
 SRCS =	CallFunction.cc				\
-	CarpetParamCheck.cc			\
 	CarpetStartup.cc			\
 	Checksum.cc				\
 	Comm.cc					\