1. Added interface for AMR. Carpet still contains no clusterer, but

1. Added interface for AMR.  Carpet still contains no clusterer, but
the application can now determine where to refine.

2. Made scalars (again) available all the time.  This allows the
scheduler to use scalars.

darcs-hash:20010826115851-07bb3-eca06f09fe49430522d9770bdf5ff2e99f7578a6.gz

5 files changed, 136 insertions, 44 deletions

diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 3357a4297..8d3d590b9 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -12,7 +12,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.4 2001/08/17 19:24:03 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.5 2001/08/26 13:58:51 schnetter Exp $";
 
 
 
@@ -22,8 +22,24 @@ namespace Carpet {
   
   
   
-  static void Recompose_gh (cGH* cgh, gh<dim>* hh);
-  static void Recompose_ga (cGH* cgh, gh<dim>* hh);
+  static bool do_recompose = false;
+  static gh<dim>::rexts next_bbsss;
+  static gh<dim>::rprocs next_pss;
+  
+  
+  
+  static void Adapt (cGH* cgh, int reflevels, gh<dim>* hh);
+  
+  
+  
+  void RegisterRecomposeRegions (const gh<dim>::rexts& bbsss,
+				 const gh<dim>::rprocs& pss)
+  {
+    // save the region information for the next regridding
+    next_bbsss = bbsss;
+    next_pss = pss;
+    do_recompose = true;
+  }
   
   
   
@@ -34,17 +50,47 @@ namespace Carpet {
     assert (component == -1);
     Checkpoint ("%*sRecompose", 2*reflevel, "");
     
-    Recompose_gh (cgh, hh);
+    // Check whether to recompose
+    if (!do_recompose) return;
     
-    // Recompose only grid functions, and not grid arrays
-    Recompose_ga (cgh, hh0);
+    // Recompose
+    hh->recompose (next_bbsss, next_pss);
     
+    if (verbose && CCTK_MyProc(cgh)==0) {
+      cout << endl;
+      cout << "New bounding boxes:" << 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:" << 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;
+    }
+    
+    // Don't recompose to these regions any more
+    do_recompose = false;
+    
+    // Adapt grid scalars
+    Adapt (cgh, hh->reflevels(), hh0);
+    
+    // Adapt grid arrays
     for (int group=0; group<CCTK_NumGroups(); ++group) {
       switch (CCTK_GroupTypeI(group)) {
       case CCTK_SCALAR:
 	break;
       case CCTK_ARRAY:
-	Recompose_ga (cgh, arrdata[group].hh);
+	Adapt (cgh, hh->reflevels(), arrdata[group].hh);
       case CCTK_GF:
 	break;
       default:
@@ -55,13 +101,12 @@ namespace Carpet {
   
   
   
-  static void Recompose_gh (cGH* cgh, gh<dim>* hh)
+  static void Adapt (cGH* cgh, const int reflevels, gh<dim>* hh)
   {
     DECLARE_CCTK_PARAMETERS;
     
     const int nprocs    = CCTK_nProcs(cgh);
-    const int reflevels = max_refinement_levels; // arbitrary value
-    const int mglevels  = 1;	// arbitrary value
+    const int mglevels  = 1;	// for now
     vector<vector<bbox<int,dim> > > bbss(reflevels);
     // note: what this routine calls "ub" is "ub+str" elsewhere
     vect<int,dim> rstr = hh->baseextent.stride();
@@ -72,20 +117,16 @@ namespace Carpet {
 	// save old values
 	const vect<int,dim> oldrlb = rlb;
 	const vect<int,dim> oldrub = rub;
-	const vect<int,dim> oldrstr = rstr;
-	// calculate extent and centre
+	// calculate extent
 	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;
+ 	// refine 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));
@@ -122,24 +163,25 @@ namespace Carpet {
   
   
   
-  static void Recompose_ga (cGH* cgh, gh<dim>* hh)
+  void MakeRegions_RefineCentre (cGH* cgh, const int reflevels,
+				 gh<dim>::rexts& bbsss, gh<dim>::rprocs& pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
     const int nprocs    = CCTK_nProcs(cgh);
-    const int reflevels = max_refinement_levels; // arbitrary value
     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;
-	const vect<int,dim> oldrstr = rstr;
 	// calculate extent and centre
 	const vect<int,dim> rextent = rub - rlb;
 	const vect<int,dim> rcentre = rlb + (rextent / 2 / rstr) * rstr;
@@ -149,10 +191,10 @@ namespace Carpet {
 	// 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;
+	// 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));
@@ -173,9 +215,9 @@ namespace Carpet {
       }
       bbss[rl] = bbs;
     }
-    vector<vector<vector<bbox<int,dim> > > > bbsss
-      = hh->make_multigrid_boxes(bbss, mglevels);
-    vector<vector<int> > pss(bbss.size());
+    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
+    
+    pss.resize(bbss.size());
     for (int rl=0; rl<reflevels; ++rl) {
       pss[rl] = vector<int>(bbss[rl].size());
       // make sure all processors have the same number of components
@@ -184,7 +226,6 @@ namespace Carpet {
 	pss[rl][c] = c % nprocs; // distribute among processors
       }
     }
-    hh->recompose(bbsss, pss);
   }
   
 } // namespace Carpet
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index a9000952b..7ae52773e 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -10,7 +10,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.5 2001/08/17 19:15:59 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.6 2001/08/26 13:59:00 schnetter Exp $";
 
 
 
@@ -227,6 +227,12 @@ namespace Carpet {
     mglevel   = 0;
     component = -1;
     
+    // Invent a refinement structure
+    gh<dim>::rexts bbsss;
+    gh<dim>::rprocs pss;
+    MakeRegions_RefineCentre (cgh, maxreflevels, bbsss, pss);
+    RegisterRecomposeRegions (bbsss, pss);
+    
     // Recompose grid hierarchy
     Recompose (cgh); 
     
@@ -242,9 +248,9 @@ namespace Carpet {
     set_component (cgh, -1);
     
     // Enable storage for all groups if desired
-    // XXX
-    if (true || enable_all_storage) {
-      for (int group=0; group<CCTK_NumGroups(); ++group) {
+    for (int group=0; group<CCTK_NumGroups(); ++group) {
+      // XXX
+      if (true || enable_all_storage || CCTK_GroupTypeI(group)==CCTK_SCALAR) {
 	EnableGroupStorage (cgh, CCTK_GroupName(group));
       }
     }
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index 5c12b3ce3..7b75869e1 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.1 2001/07/09 09:00:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.2 2001/08/26 13:59:08 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
@@ -88,6 +88,8 @@ namespace Carpet {
 #include "carpet_public.h"
   }
   
+  
+  
   // Registered functions
   void* SetupGH (tFleshConfig* fc, int convLevel, cGH* cgh);
   
@@ -113,6 +115,15 @@ namespace Carpet {
   
   
   
+  // Functions for recomposing the grid hierarchy
+  void RegisterRecomposeRegions (const gh<dim>::rexts& bbsss,
+				 const gh<dim>::rprocs& pss);
+  
+  void MakeRegions_RefineCentre (cGH* cgh, int reflevels,
+				 gh<dim>::rexts& bbsss, gh<dim>::rprocs& pss);
+  
+  
+  
   // Helper functions
   void set_reflevel (cGH* cgh, int rl);
   void set_mglevel (cGH* cgh, int ml);
@@ -162,12 +173,13 @@ namespace Carpet {
   
   // Component iterator
   
-#define BEGIN_COMPONENT_LOOP(cgh)		\
-  do {						\
-    int _cl;					\
-    assert (component==-1);			\
-    set_component ((cgh), 0);			\
-    for (;;) {					\
+#define BEGIN_COMPONENT_LOOP(cgh)			\
+  do {							\
+    int _cl;						\
+    assert (reflevel>=0 && reflevel<hh->reflevels());	\
+    assert (component==-1);				\
+    set_component ((cgh), 0);				\
+    for (;;) {						\
       {
 #define END_COMPONENT_LOOP(cgh)				\
       }							\
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index c19ec197f..5207ebd62 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.4 2001/08/23 17:52:58 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.5 2001/08/26 13:59:15 schnetter Exp $";
 
 
 
@@ -182,9 +182,12 @@ namespace Carpet {
   
   void set_component (cGH* cgh, const int c)
   {
+    assert (reflevel>=0 && reflevel<hh->reflevels());
     assert (c==-1 || (c>=0 && c<hh->components(reflevel)));
     component = c;
-      
+    assert (component==-1
+	    || (mglevel>=0 && mglevel<hh->mglevels(reflevel,component)));
+    
     if (component == -1) {
       // Global mode -- no component is active
       
@@ -218,8 +221,35 @@ namespace Carpet {
 	  const int var   = n - CCTK_FirstVarIndexI(group);
 	  assert (var>=0);
 	  
-	  // Scalars, arrays, and grid functions cannot be accessed
-	  cgh->data[n][tl] = 0;
+	  if (CCTK_GroupTypeI(group) != CCTK_SCALAR) {
+	    // Arrays and grid functions cannot be accessed
+	    
+	    cgh->data[n][tl] = 0;
+	    
+	  } else {
+	    // Scalars can be accessed
+	    
+	    if (CCTK_QueryGroupStorageI(cgh, group)) {
+	      // Group has storage
+	      
+	      assert (group<(int)arrdata.size());
+	      assert (var<(int)arrdata[group].data.size());
+	      assert (arrdata[group].data[var]);
+	      const int c = CCTK_MyProc(cgh);
+	      assert (hh->is_local(reflevel,c));
+	      cgh->data[n][tl]
+		= ((*arrdata[group].data[var])
+		   (-tl, reflevel, c, mglevel)->storage());
+	      assert (cgh->data[n][tl]);
+	      
+	    } else {
+	      // Group has no storage
+	      
+	      cgh->data[n][tl] = 0;
+	      
+	    } // if ! has storage
+	    
+	  } // if group type is SCALAR
 	  
 	} // for tl
       }	// for n
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index ac868e5a9..f862195d3 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.7 2001/07/04 12:29:51 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.8 2001/08/26 13:59:31 schnetter Exp $
 
  ***************************************************************************/
 
@@ -79,7 +79,10 @@ ostream& output (ostream& os, const vector<T>& v) {
 template ostream& output (ostream& os, const list<bbox<int,3> >& l);
 template ostream& output (ostream& os, const set<bbox<int,3> >& s);
 template ostream& output (ostream& os, const set<bboxset<int,3> >& s);
+template ostream& output (ostream& os, const vector<int>& v);
+template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
 template ostream& output (ostream& os, const vector<list<bbox<int,3> > >& v);
+template ostream& output (ostream& os, const vector<vector<int> >& v);
 template ostream& output (ostream& os, const vector<vector<bbox<int,3> > >& v);
 template ostream& output (ostream& os, const vector<vector<vector<bbox<int,3> > > >& v);
 #endif