Two refinement levels now work, but only without symmetry boundaries.

Two refinement levels now work, but only without symmetry boundaries.
"work" is here defined by the eyeball.  With symmetry boundaries,
there is either a bug or a bad interaction of restricting,
prolongating, and the symmetry boundary on the coarse grid.

darcs-hash:20010313164030-f6438-2c18a10421487f6328e0e782322332c67fb203fd.gz

6 files changed, 68 insertions, 39 deletions

diff --git a/Carpet/Carpet/src/carpet.cc b/Carpet/Carpet/src/carpet.cc
index 3ad87c574..1ebb4a70d 100644
--- a/Carpet/Carpet/src/carpet.cc
+++ b/Carpet/Carpet/src/carpet.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.8 2001/03/13 13:06:52 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.9 2001/03/13 17:40:30 eschnett 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,7 +32,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.8 2001/03/13 13:06:52 eschnett Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.9 2001/03/13 17:40:30 eschnett Exp $";
 
 
 
@@ -241,13 +241,13 @@ namespace Carpet {
     
     // initialise cgh
     cgh->cctk_convlevel = mglevel;
-    for (int d=0; d<dim; ++d) {
-      cgh->cctk_levfac[d] = 1;
-    }
     
     // Recompose grid hierarchy
     Recompose (cgh); 
     
+    // enace current position
+    enact_reflevel (cgh);
+    
     Checkpoint ("done with SetupGH.");
     
     // We register only once, ergo we get only one handle, ergo there
@@ -507,6 +507,7 @@ namespace Carpet {
   
   
   
+#if 0
   int ScheduleTraverse (cGH* cgh, const char* rfrPoint)
   {
     // traverse all functions on all refinement levels on all
@@ -526,6 +527,7 @@ namespace Carpet {
     
     return 0;
   }
+#endif
   
   
   
@@ -554,7 +556,7 @@ namespace Carpet {
       // set Cactus parameters to pseudo values
       for (int d=0; d<dim; ++d) {
 	cgh->cctk_lsh[d]      = 0xdeadbeef;
-	cgh->cctk_gsh[d]      = 0xdeadbeef;
+// 	cgh->cctk_gsh[d]      = 0xdeadbeef;
 	cgh->cctk_bbox[2*d  ] = 0xdeadbeef;
 	cgh->cctk_bbox[2*d+1] = 0xdeadbeef;
 	cgh->cctk_lbnd[d]     = 0xdeadbeef;
@@ -643,21 +645,26 @@ namespace Carpet {
 	  
 	  // set Cactus parameters
 	  for (int d=0; d<dim; ++d) {
-	    typedef bbox<int,dim> ibbox;
-	    const ibbox ext = dd->boxes[reflevel][component][mglevel].exterior;
-	    const ibbox base = hh->baseextent;
+	    const bbox<int,dim>& ext
+	      = dd->boxes[reflevel][component][mglevel].exterior;
+	    const bbox<int,dim>& base = hh->baseextent;
 	    cgh->cctk_lsh[d] = (ext.shape() / ext.stride())[d];
-	    cgh->cctk_gsh[d]
-	      = ((base.shape() / base.stride() + dd->lghosts + dd->ughosts)[d]
-		 * cgh->cctk_levfac[d]);
+// 	    cgh->cctk_gsh[d]
+// 	      = ((base.shape() / base.stride() + dd->lghosts + dd->ughosts)[d]
+// 		 * cgh->cctk_levfac[d]);
 	    cgh->cctk_lbnd[d] = (ext.lower() / ext.stride())[d];
 	    cgh->cctk_ubnd[d] = (ext.upper() / ext.stride())[d];
-// 	    cgh->cctk_bbox[2*d  ]
-// 	      = reflevel==0 && cgh->cctk_lbnd[d] == 0;
-// 	    cgh->cctk_bbox[2*d+1]
-// 	      = reflevel==0 && cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1;
-	    cgh->cctk_bbox[2*d  ] = (ext.lower() < base.lower())[d];
-	    cgh->cctk_bbox[2*d+1] = (ext.upper() > base.upper())[d];
+	    // no outer boundaries on the finer grids
+	    cgh->cctk_bbox[2*d  ]
+	      = reflevel==0 && cgh->cctk_lbnd[d] == 0;
+	    cgh->cctk_bbox[2*d+1]
+	      = reflevel==0 && cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1;
+#if 0
+	    // do allow outer boundaries on the finer grids
+	    // (but this is generally inconsistent -- c. f. periodicity)
+ 	    cgh->cctk_bbox[2*d  ] = (ext.lower() < base.lower())[d];
+ 	    cgh->cctk_bbox[2*d+1] = (ext.upper() > base.upper())[d];
+#endif
 	    for (int stg=0; stg<CCTK_NSTAGGER; ++stg) {
 	      // TODO: support staggering
 	      cgh->cctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d];
@@ -1049,12 +1056,25 @@ namespace Carpet {
     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 + rstr;
+	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) the center only
-	rlb /= hh->reffact;
-	rub /= hh->reffact;
+	// refine (arbitrarily) around the center only
+	rlb = rcentre - (newrextent/2 / rstr) * rstr;
+	rub = rlb + newrextent - rstr;
+// 	// refine (arbitrarily) around the lower boundary only
+// 	rlb = rlb;
+// 	rub = rlb + newrextent - rstr;
+	assert (all(rlb >= oldrlb && rub <= oldrub));
       }
       vector<bbox<int,dim> > bbs(nprocs);
       for (int c=0; c<nprocs; ++c) {
@@ -1340,9 +1360,13 @@ namespace Carpet {
     assert (component == -1);
     
     // Change
+    const bbox<int,dim>& base = hh->baseextent;
     reflevelfactor = (int)floor(pow(hh->reffact, reflevel)+0.5);
     cgh->cctk_delta_time = base_delta_time / reflevelfactor;
     for (int d=0; d<dim; ++d) {
+      cgh->cctk_gsh[d]
+	= ((base.shape() / base.stride() + dd->lghosts + dd->ughosts)[d]
+	   * reflevelfactor);
       cgh->cctk_levfac[d] = reflevelfactor;
     }
   }
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index 31405e865..5b4d10bfc 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.5 2001/03/13 13:06:59 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.6 2001/03/13 17:40:37 eschnett Exp $
 
 #include <cassert>
 #include <cstdio>
@@ -443,17 +443,17 @@ template<int outdim>
 int IOASCII<outdim>::GetGridOffset (cGH* cgh, int dir,
 				    const char* itempl, const char* iglobal,
 				    const char* ctempl, const char* cglobal,
-				    double cfallback)
+				    CCTK_REAL cfallback)
 {
   /* First choice: explicit coordinate */
   char cparam[1000];
   sprintf (cparam, ctempl, outdim);
   if (CCTK_ParameterQueryTimesSet (cparam, CCTK_THORNSTRING) > 0) {
     int ptype;
-    const double* const pcoord
-      = (double*)CCTK_ParameterGet (cparam, CCTK_THORNSTRING, &ptype);
+    const CCTK_REAL* const pcoord
+      = (CCTK_REAL*)CCTK_ParameterGet (cparam, CCTK_THORNSTRING, &ptype);
     assert (pcoord);
-    const double coord = *pcoord;
+    const CCTK_REAL coord = *pcoord;
     assert (ptype == PARAMETER_REAL);
     return CoordToOffset (cgh, dir, coord);
   }
@@ -474,10 +474,10 @@ int IOASCII<outdim>::GetGridOffset (cGH* cgh, int dir,
   /* Third choice: explicit global coordinate */
   if (CCTK_ParameterQueryTimesSet (cglobal, "IO") > 0) {
     int ptype;
-    const double* const pcoord
-      = (double*)CCTK_ParameterGet (cglobal, "IO", &ptype);
+    const CCTK_REAL* const pcoord
+      = (CCTK_REAL*)CCTK_ParameterGet (cglobal, "IO", &ptype);
     assert (pcoord);
-    const double coord = *pcoord;
+    const CCTK_REAL coord = *pcoord;
     assert (ptype == PARAMETER_REAL);
     return CoordToOffset (cgh, dir, coord);
   }
@@ -500,7 +500,7 @@ int IOASCII<outdim>::GetGridOffset (cGH* cgh, int dir,
 
 
 template<int outdim>
-int IOASCII<outdim>::CoordToOffset (cGH* cgh, int dir, double coord)
+int IOASCII<outdim>::CoordToOffset (cGH* cgh, int dir, CCTK_REAL coord)
 {
   CCTK_REAL lower, upper;
   CCTK_CoordRange (cgh, &lower, &upper, dir, 0, "cart3d");
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index c6b4ad3d6..ba3cb79df 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.6 2001/03/13 13:06:59 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.7 2001/03/13 17:40:38 eschnett Exp $
 
  ***************************************************************************/
 
@@ -30,7 +30,7 @@
 #  include "dh.hh"
 #endif
 
-#define DEBUG_OUTPUT
+#undef DEBUG_OUTPUT
 
 
 
diff --git a/CarpetExtra/IDScalarWave/schedule.ccl b/CarpetExtra/IDScalarWave/schedule.ccl
index cee342d08..d4736d90b 100644
--- a/CarpetExtra/IDScalarWave/schedule.ccl
+++ b/CarpetExtra/IDScalarWave/schedule.ccl
@@ -1,12 +1,12 @@
 # Schedule definitions for thorn IDScalarWave
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWave/schedule.ccl,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWave/schedule.ccl,v 1.2 2001/03/13 17:40:38 eschnett Exp $
 
-schedule IDScalarWave_CheckParameters at CCTK_PARAMCHECK
+schedule IDScalarWave_CheckParameters at PARAMCHECK
 {
   LANG: Fortran
 } "Check parameters"
 
-schedule IDScalarWave_InitialData at CCTK_INITIAL
+schedule IDScalarWave_InitialData as WaveToy_InitialData at INITIAL
 {
   STORAGE:       wavetoy::scalarevolve
   LANG:          Fortran
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_rl2.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_rl2.par
index d53c7890f..9b28340a4 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_rl2.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_rl2.par
@@ -7,11 +7,11 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/Attic/wavetoyf77_rad_rl2.par,v 1.2 2001/03/13 13:07:32 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/Attic/wavetoyf77_rad_rl2.par,v 1.3 2001/03/13 17:40:38 eschnett Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetSlab IDScalarWave WaveToyF77"
 
-Cactus::cctk_itlast		= 1 # 120
+Cactus::cctk_itlast		= 120
  
 Time::dtfac			= 0.5
 
@@ -24,7 +24,7 @@ Carpet::max_refinement_levels	= 2
 grid::type			= byspacing
 grid::mode			= octant
 grid::dxyz			= 0.3
- 
+
 IO::outdir			= "wavetoyf77_rad_rl2"
 
 IOBasic::outinfo_every		= 1 # 10
diff --git a/CarpetExtra/WaveToyF77/schedule.ccl b/CarpetExtra/WaveToyF77/schedule.ccl
index 70a324fb7..fdb922a8c 100644
--- a/CarpetExtra/WaveToyF77/schedule.ccl
+++ b/CarpetExtra/WaveToyF77/schedule.ccl
@@ -1,5 +1,5 @@
 # Schedule definitions for thorn WaveToy77
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/schedule.ccl,v 1.2 2001/03/07 13:01:47 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/schedule.ccl,v 1.3 2001/03/13 17:40:38 eschnett Exp $
 
 STORAGE: scalarevolve
 
@@ -20,6 +20,11 @@ schedule WaveToyF77_Evolution as WaveToy_Evolution at EVOL
   SYNC: scalarevolve
 } "Evolution of 3D wave equation"
 
+schedule WaveToyF77_Boundaries as WaveToy_Boundaries at INITIAL after WaveToy_InitialData
+{
+  LANG: Fortran
+} "Boundaries of 3D wave equation"
+
 schedule WaveToyF77_Boundaries as WaveToy_Boundaries at EVOL after WaveToy_Evolution
 {
   LANG: Fortran