CarpetRegrid2: Correct checkpointing error. Implement cell centering.

Correct error with checkpointing and recovery.  Store the last
regridding iteration in grid variables, so that it is conserved during
checkpointing and recovery.

Implement cell centered refinement.  Introduce a new parameter
"snap_to_coarse" that makes the fine grid always cover complete coarse
grid cells.

Add a parameter "min_fraction" which specifies when to combine several
small grids to form a single, larger, combined grid.

darcs-hash:20070112233156-dae7b-3f32c9558e29f25f9bc1ba533985a4392877feab.gz

4 files changed, 82 insertions, 17 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 95f551a69..206c7c7bc 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -5,6 +5,7 @@ IMPLEMENTS: CarpetRegrid2
 USES INCLUDE HEADER: bbox.hh
 USES INCLUDE HEADER: bboxset.hh
 USES INCLUDE HEADER: defs.hh
+USES INCLUDE HEADER: dh.hh
 USES INCLUDE HEADER: gh.hh
 USES INCLUDE HEADER: region.hh
 USES INCLUDE HEADER: vect.hh
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index df17204dd..dbcdee93b 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -6,27 +6,24 @@ BOOLEAN verbose "Display regridding information on the terminal" STEERABLE=alway
 
 
 
-CCTK_INT num_centres "Number of refinement centres" STEERABLE=always
-{
-  0:3 :: ""
-} 0
-
-
-
 CCTK_INT min_distance "Minimum distance (in grid points) between coarse and fine grid boundaries" STEERABLE=always
 {
   0:* :: ""
 } 4
 
-
-
-CCTK_REAL min_fraction "Minimum fraction of required refined points that need to be present in a refined region"
+CCTK_REAL min_fraction "Minimum fraction of required refined points that need to be present in a refined region" STEERABLE=always
 {
   0:* :: ""
 } 0.9
 
 
 
+BOOLEAN snap_to_coarse "Ensure that the fine grid extent coincides with coarse grid points" STEERABLE=always
+{
+} "no"
+
+
+
 CCTK_INT regrid_every "Regrid every n time steps" STEERABLE=always
 {
   -1   :: "regrid never"
@@ -46,6 +43,13 @@ BOOLEAN symmetry_rotating180 "Ensure a 180 degree rotating symmetry about the z
 
 
 
+CCTK_INT num_centres "Number of refinement centres" STEERABLE=always
+{
+  0:3 :: ""
+} 0
+
+
+
 ################################################################################
 
 CCTK_INT num_levels_1 "Number of refinement levels for this centre" STEERABLE=always
diff --git a/Carpet/CarpetRegrid2/schedule.ccl b/Carpet/CarpetRegrid2/schedule.ccl
index bdc92e106..6299ba209 100644
--- a/Carpet/CarpetRegrid2/schedule.ccl
+++ b/Carpet/CarpetRegrid2/schedule.ccl
@@ -10,7 +10,7 @@ SCHEDULE CarpetRegrid2_ParamCheck AT paramcheck
   LANG: C
 } "Check parameters"
 
-SCHEDULE CarpetRegrid2_Initialise AT initial
+SCHEDULE CarpetRegrid2_Initialise AT wragh
 {
   LANG: C
   OPTIONS: global
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 96b1e16a1..1feb322c4 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -74,13 +74,43 @@ namespace CarpetRegrid2 {
              rvect const & origin, rvect const & scale,
              gh const & hh, int const rl)
   {
+    DECLARE_CCTK_PARAMETERS;
+    
     ivect const levfac = hh.reffacts.at(rl);
     assert (all (hh.baseextent.stride() % levfac == 0));
     ivect const istride = hh.baseextent.stride() / levfac;
     
-    return (ivect (floor ((rpos - origin) * scale / rvect(istride) +
-                          static_cast<CCTK_REAL> (0.5))) *
-            istride);
+    ivect ipos = (ivect (floor ((rpos - origin) * scale / rvect(istride) +
+                                static_cast<CCTK_REAL> (0.5))) *
+                  istride);
+    
+    if (snap_to_coarse) {
+      if (rl > 0) {
+        ivect const clevfac = hh.reffacts.at(rl - 1);
+        assert (all (hh.baseextent.stride() % clevfac == 0));
+        ivect const cistride = hh.baseextent.stride() / clevfac;
+        ipos = ipos / cistride * cistride;
+      }
+    }
+    
+    if (hh.refcent == cell_centered) {
+#if 0
+      if (rl > 0) {
+        ivect const clevfac = hh.reffacts.at(rl - 1);
+        assert (all (hh.baseextent.stride() % clevfac == 0));
+        ivect const cistride = hh.baseextent.stride() / clevfac;
+        
+        ivect const offset = (cistride / istride + 1) % 2;
+        
+        assert (all (istride % offset == 0));
+        ipos += offset;
+      }
+#endif
+      assert (all (istride % 2 == 0));
+      ipos += istride / 2;
+    }
+    
+    return ipos;
   }
   
   // Convert a coordinate location to an index location, rounding in
@@ -90,13 +120,43 @@ namespace CarpetRegrid2 {
               rvect const & origin, rvect const & scale,
               gh const & hh, int const rl)
   {
+    DECLARE_CCTK_PARAMETERS;
+    
     ivect const levfac = hh.reffacts.at(rl);
     assert (all (hh.baseextent.stride() % levfac == 0));
     ivect const istride = hh.baseextent.stride() / levfac;
     
-    return (ivect (ceil ((rpos - origin) * scale / rvect(istride) -
-                         static_cast<CCTK_REAL> (0.5))) *
-            istride);
+    ivect ipos = (ivect (ceil ((rpos - origin) * scale / rvect(istride) -
+                               static_cast<CCTK_REAL> (0.5))) *
+                  istride);
+    
+    if (snap_to_coarse) {
+      if (rl > 0) {
+        ivect const clevfac = hh.reffacts.at(rl - 1);
+        assert (all (hh.baseextent.stride() % clevfac == 0));
+        ivect const cistride = hh.baseextent.stride() / clevfac;
+        ipos = (ipos + cistride - 1) / cistride * cistride;
+      }
+    }
+    
+    if (hh.refcent == cell_centered) {
+#if 0
+      if (rl > 0) {
+        ivect const clevfac = hh.reffacts.at(rl - 1);
+        assert (all (hh.baseextent.stride() % clevfac == 0));
+        ivect const cistride = hh.baseextent.stride() / clevfac;
+        
+        ivect const offset = (cistride / istride + 1) % 2;
+        
+        assert (all (istride % offset == 0));
+        ipos -= offset;
+      }
+#endif
+      assert (all (istride % 2 == 0));
+      ipos -= istride / 2;
+    }
+    
+    return ipos;
   }
   
   // Convert an index location to a coordinate location