Added box-in-box capability to Carpet.

darcs-hash:20011205023155-07bb3-2b8f43070aac0e38ad0f896b6ae73c8baa7608af.gz

8 files changed, 370 insertions, 21 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index 176b0e182..be9b780ea 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn Carpet
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.10 2001/07/11 17:41:13 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.11 2001/12/05 03:31:55 schnetter Exp $
 
 shares: Cactus
 
@@ -125,3 +125,110 @@ BOOLEAN veryverbose "Display a lot of info on the screen"
 BOOLEAN barriers "Insert barriers at strategic places for debugging purposes (slows down execution)"
 {
 } "no"
+
+
+
+KEYWORD processor_topology "How to determine the processor topology"
+{
+  "manual"    :: "Specified by processor_topology_*"
+  "automatic" :: "Choose the topology automatically"
+} "automatic"
+
+CCTK_INT processor_topology_3d_x "Number of processors in x-direction"
+{
+  1:* :: "must be positive"
+} 1
+
+CCTK_INT processor_topology_3d_y "Number of processors in y-direction"
+{
+  1:* :: "must be positive"
+} 1
+
+CCTK_INT processor_topology_3d_z "Number of processors in z-direction"
+{
+  1:* :: "must be positive"
+} 1
+
+
+
+CCTK_INT l1xmin "Lower boundary of level 1 box in x-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l1ymin "Lower boundary of level 1 box in y-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l1zmin "Lower boundary of level 1 box in z-direction"
+{
+  : :: ""
+} 0
+
+CCTK_INT l1xmax "Upper boundary of level 1 box in x-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l1ymax "Upper boundary of level 1 box in y-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l1zmax "Upper boundary of level 1 box in z-direction"
+{
+  : :: ""
+} -1
+
+
+
+CCTK_INT l2xmin "Lower boundary of level 2 box in x-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l2ymin "Lower boundary of level 2 box in y-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l2zmin "Lower boundary of level 2 box in z-direction"
+{
+  : :: ""
+} 0
+
+CCTK_INT l2xmax "Upper boundary of level 2 box in x-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l2ymax "Upper boundary of level 2 box in y-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l2zmax "Upper boundary of level 2 box in z-direction"
+{
+  : :: ""
+} -1
+
+
+
+CCTK_INT l3xmin "Lower boundary of level 3 box in x-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l3ymin "Lower boundary of level 3 box in y-direction"
+{
+  : :: ""
+} 0
+CCTK_INT l3zmin "Lower boundary of level 3 box in z-direction"
+{
+  : :: ""
+} 0
+
+CCTK_INT l3xmax "Upper boundary of level 3 box in x-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l3ymax "Upper boundary of level 3 box in y-direction"
+{
+  : :: ""
+} -1
+CCTK_INT l3zmax "Upper boundary of level 3 box in z-direction"
+{
+  : :: ""
+} -1
diff --git a/Carpet/Carpet/schedule.ccl b/Carpet/Carpet/schedule.ccl
index 660ea38fa..0a4778ea0 100644
--- a/Carpet/Carpet/schedule.ccl
+++ b/Carpet/Carpet/schedule.ccl
@@ -1,5 +1,10 @@
 # Schedule definitions for thorn Carpet
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/schedule.ccl,v 1.2 2001/03/07 12:59:40 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/schedule.ccl,v 1.3 2001/12/05 03:31:55 schnetter Exp $
+
+schedule CarpetParamCheck at PARAMCHECK
+{
+  LANG: C
+} "Check parameters"
 
 schedule CarpetStartup at STARTUP as Driver_Startup
 {
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 557b77377..72b80c51d 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -1,4 +1,5 @@
 #include <assert.h>
+#include <math.h>
 #include <stdlib.h>
 
 #include <vector>
@@ -12,7 +13,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.8 2001/11/14 17:57:55 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.9 2001/12/05 03:31:56 schnetter Exp $";
 
 
 
@@ -156,12 +157,11 @@ namespace Carpet {
   // 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, gh<dim>::rprocs& pss)
+				 gh<dim>::rexts& bbsss)
   {
     DECLARE_CCTK_PARAMETERS;
     
-    const int nprocs    = CCTK_nProcs(cgh);
-    const int mglevels  = 1;	// 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
@@ -191,6 +191,81 @@ namespace Carpet {
 	// require rub<oldrub because we really want rub-rstr<=oldrub-oldstr
 	assert (all(rlb >= oldrlb && rub < oldrub));
       }
+      vector<bbox<int,dim> > bbs(1);
+      bbs[0] = bbox<int,dim>(rlb, rub-rstr, rstr);
+      bbss[rl] = bbs;
+    }
+    bbsss = hh->make_multigrid_boxes(bbss, mglevels);
+  }
+  
+  
+  
+  void MakeRegions_AsSpecified (const cGH* cgh, const int reflevels,
+				gh<dim>::rexts& bbsss)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    const int mglevels = 1;	// arbitrary value
+    
+    // note: what this routine calls "ub" is "ub+str" elsewhere
+    const vect<int,dim> rstr = hh->baseextent.stride();
+    const vect<int,dim> rlb  = hh->baseextent.lower();
+    const vect<int,dim> rub  = hh->baseextent.upper() + rstr;
+    
+    if (reflevels>3) {
+      CCTK_WARN (0, "Cannot currently specify refinement regions for more than 3 refinement levels");
+    }
+    
+    assert (reflevels<4);
+    vector<vect<int,dim> > lower(4), upper(4);
+    lower[0] = rlb;
+    upper[0] = rub;
+    lower[1] = rstr * vect<int,dim> (l1xmin, l1ymin, l1zmin);
+    upper[1] = rstr * vect<int,dim> (l1xmax, l1ymax, l1zmax);
+    lower[2] = rstr * vect<int,dim> (l2xmin, l2ymin, l2zmin);
+    upper[2] = rstr * vect<int,dim> (l2xmax, l2ymax, l2zmax);
+    lower[3] = rstr * vect<int,dim> (l3xmin, l3ymin, l3zmin);
+    upper[3] = rstr * vect<int,dim> (l3xmax, l3ymax, l3zmax);
+    
+    vector<vector<bbox<int,dim> > > bbss(reflevels);
+    
+    for (int rl=0; rl<reflevels; ++rl) {
+      const int fact = floor(pow(hh->reffact, rl) + 0.5);
+      assert (all (rstr % fact == 0));
+      const vect<int,dim> str (rstr / fact);
+      const vect<int,dim> lb  (lower[rl]);
+      const vect<int,dim> ub  (upper[rl]);
+      assert (all(lb<=ub && lb>=rlb && ub-str<=rub-rstr
+		  && lb%str==0 && ub%str==0));
+      vector<bbox<int,dim> > bbs(1);
+      bbs[0] = bbox<int,dim>(lb, ub-str, str);
+      bbss[rl] = bbs;
+    }
+    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)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    const int nprocs    = CCTK_nProcs(cgh);
+    const int mglevels  = 1;	// arbitrary value
+    
+    vector<vector<bbox<int,dim> > > bbss(bbsss.size());
+    
+    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;
@@ -210,13 +285,83 @@ namespace Carpet {
       bbss[rl] = bbs;
     }
     bbsss = hh->make_multigrid_boxes(bbss, mglevels);
+  }
+  
+  
+  
+  void SplitRegions_AsSpecified (const cGH* cgh, gh<dim>::rexts& bbsss)
+  {
+    DECLARE_CCTK_PARAMETERS;
     
-    pss.resize(bbss.size());
-    for (int rl=0; rl<reflevels; ++rl) {
-      pss[rl] = vector<int>(bbss[rl].size());
+    const int nprocs    = CCTK_nProcs(cgh);
+    const int mglevels  = 1;	// arbitrary value
+    
+    vector<vector<bbox<int,dim> > > bbss(bbsss.size());
+    
+    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);
+	  }
+	}
+      }
+      bbss[rl] = bbs;
+    }
+    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 MakeProcessors_RoundRobin (const cGH* cgh, const gh<dim>::rexts& bbsss,
+				  gh<dim>::rprocs& pss)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    const int nprocs    = CCTK_nProcs(cgh);
+    
+    pss.resize(bbsss.size());
+    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
+      pss[rl] = vector<int>(bbsss[rl].size());
       // make sure all processors have the same number of components
-      assert (bbss[rl].size() % nprocs == 0);
-      for (int c=0; c<(int)bbss[rl].size(); ++c) {
+      assert (bbsss[rl].size() % nprocs == 0);
+      for (int c=0; c<(int)bbsss[rl].size(); ++c) {
 	pss[rl][c] = c % nprocs; // distribute among processors
       }
     }
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index d8e688c3e..2ff02193f 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.11 2001/11/15 16:41:32 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.12 2001/12/05 03:31:56 schnetter Exp $";
 
 
 
@@ -38,8 +38,7 @@ namespace Carpet {
     
     // Refinement information
     maxreflevels = max_refinement_levels;
-    maxreflevelfact
-      = (int)floor(pow((double)refinement_factor, maxreflevels-1) + 0.5);
+    maxreflevelfact = floor(pow(refinement_factor, maxreflevels-1) + 0.5);
     
     // Ghost zones
     vect<int,dim> lghosts, ughosts;
@@ -235,7 +234,16 @@ namespace Carpet {
     // Invent a refinement structure
     gh<dim>::rexts bbsss;
     gh<dim>::rprocs pss;
-    MakeRegions_RefineCentre (cgh, maxreflevels, bbsss, pss);
+//     MakeRegions_RefineCentre (cgh, maxreflevels, bbsss);
+    MakeRegions_AsSpecified (cgh, maxreflevels, 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);
     RegisterRecomposeRegions (bbsss, pss);
     
     // Recompose grid hierarchy
diff --git a/Carpet/Carpet/src/carpet_public.h b/Carpet/Carpet/src/carpet_public.h
index 168ef725b..c80d4d3d8 100644
--- a/Carpet/Carpet/src/carpet_public.h
+++ b/Carpet/Carpet/src/carpet_public.h
@@ -1,10 +1,11 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.h,v 1.1 2001/07/09 09:00:13 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.h,v 1.2 2001/12/05 03:31:56 schnetter Exp $ */
 
 #include <mpi.h>
 
 #include "cctk.h"
 
 /* Scheduled functions */
+int CarpetParamCheck (void);
 int CarpetStartup (void);
 
 /* Helper functions */
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index 9cdc73379..9497e32aa 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.6 2001/11/05 17:53:02 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.7 2001/12/05 03:31:57 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
@@ -119,8 +119,14 @@ 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, 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);
   
   
   
diff --git a/Carpet/Carpet/src/make.code.defn b/Carpet/Carpet/src/make.code.defn
index 0c1146197..446f04cc6 100644
--- a/Carpet/Carpet/src/make.code.defn
+++ b/Carpet/Carpet/src/make.code.defn
@@ -1,8 +1,9 @@
 # Main make.code.defn file for thorn Carpet   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.2 2001/07/04 12:29:48 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.3 2001/12/05 03:31:57 schnetter Exp $
 
 # Source files in this directory
 SRCS =	CallFunction.cc				\
+	CarpetParamCheck.cc			\
 	CarpetStartup.cc			\
 	Checksum.cc				\
 	Comm.cc					\
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index f1ee65591..5905d14b9 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.4 2001/03/27 22:26:31 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.5 2001/12/05 03:31:57 schnetter Exp $
 
  ***************************************************************************/
 
@@ -70,6 +70,11 @@ public:
     elt[0]=x; elt[1]=y; elt[2]=z;
   }
   
+  vect (const T x, const T y, const T z, const T t) {
+    assert (D==4);
+    elt[0]=x; elt[1]=y; elt[2]=z; elt[3]=t;
+  }
+  
   vect (const T* const x) {
     for (int d=0; d<D; ++d) elt[d]=x[d];
   }
@@ -405,6 +410,36 @@ public:
   };
 #endif
   
+  // Higher order functions
+  vect map (T (* const func)(T x)) const {
+    vect r;
+    for (int d=0; d<D; ++d) r[d] = func(elt[d]);
+    return r;
+  }
+  
+  T fold (T (* const func)(T val, T x), T val) const {
+    for (int d=0; d<D; ++d) val = func(val, elt[d]);
+    return val;
+  }
+  
+  vect scan0 (T (* const func)(T val, T x), T val) const {
+    vect r;
+    for (int d=0; d<D; ++d) {
+      r[d] = val;
+      val = func(val, elt[d]);
+    }
+    return r;
+  }
+  
+  vect scan1 (T (* const func)(T val, T x), T val) const {
+    vect r;
+    for (int d=0; d<D; ++d) {
+      val = func(val, elt[d]);
+      r[d] = val;
+    }
+    return r;
+  }
+  
   void output (ostream& os) const;
 };
 
@@ -501,6 +536,47 @@ inline T sum (const vect<T,D>& a) {
   return r;
 }
 
+// Higher order functions
+template<class T,class TT,int D>
+inline vect<TT,D> map (TT (* const func)(T x), const vect<T,D>& a) {
+  vect<TT,D> r;
+  for (int d=0; d<D; ++d) r[d] = func(a[d]);
+  return r;
+}
+
+template<class T,class TT,int D>
+inline vect<TT,D> fold (TT (* const func)(TT val, T x), TT val,
+			const vect<T,D>& a)
+{
+  for (int d=0; d<D; ++d) val = func(val, a[d]);
+  return val;
+}
+
+template<class T,class TT,int D>
+inline vect<TT,D> scan0 (TT (* const func)(TT val, T x), TT val,
+			 const vect<T,D>& a)
+{
+  vect<TT,D> r;
+  for (int d=0; d<D; ++d) {
+    r[d] = val;
+    val = func(val, a[d]);
+  }
+  return r;
+}
+
+template<class T,class TT,int D>
+inline vect<TT,D> scan1 (TT (* const func)(TT val, T x), TT val,
+			 const vect<T,D>& a)
+{
+  vect<TT,D> r;
+  for (int d=0; d<D; ++d) {
+    val = func(val, a[d]);
+    r[d] = val;
+  }
+  return r;
+}
+
+
 
 template<class T,int D>
 inline ostream& operator<< (ostream& os, const vect<T,D>& a) {