Carpet: Use outer buffer zones

Deprecated the previous buffer zone mechanism.

Introduce new parameters Carpet::use_outer_buffer_zones and
Carpet::num_integrator_substeps.  The number of buffer zones are then
calculated automatically, and are added to the outer edge of the grid
instead of reducing the domain.

Adapt to the API change of the dh class.

Both new and old style buffer zones can exist at the same time,
although I would not do that.

darcs-hash:20050809091707-891bb-c25d68de3b04b03c06078aaee39f389b67f1416a.gz

3 files changed, 42 insertions, 12 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index 44d62df1d..3bd189e38 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -164,13 +164,24 @@ CCTK_INT prolongation_order_time "Order of prolongation operator in time" STEERA
 
 
 
-CCTK_INT buffer_width "Width of the buffer zone inside the fine grid" STEERABLE=recover
+CCTK_INT buffer_width "Width of the buffer zone inside the fine grid (DEPRECATED)" STEERABLE=recover
 {
   0:* :: "Should be the radius of the numerical domain of dependence of the time integrator, minus the number of ghost zones"
 } 0
 
 
 
+BOOLEAN use_outer_buffer_zones "Add buffer zones at the outer edge of the refined regions; the buffer width is nghostzones * num_integrator_substeps + buffer_width"
+{
+} no
+
+CCTK_INT num_integrator_substeps "Number of substeps of the time integrator"
+{
+  0:* :: ""
+} 1
+
+
+
 CCTK_STRING base_extents "Extents of base grid components, in grid point units of the finest level"
 {
   "^$" :: "leave empty for one grid component covering the whole region (default)"
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 79c7314b5..48bc71c58 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -306,7 +306,10 @@ namespace Carpet {
                       CCTK_INT convergence_factor, CCTK_INT refinement_factor,
                       CCTK_STRING base_outerbounds, CCTK_STRING base_extents,
                       CCTK_INT max_refinement_levels,
-                      CCTK_INT buffer_width, CCTK_INT prolongation_order_space,
+                      CCTK_INT prolongation_order_space,
+                      CCTK_INT buffer_width,
+                      CCTK_INT use_outer_buffer_zones,
+                      CCTK_INT num_integrator_substeps,
                       ivect & lghosts, ivect & ughosts, ivect & npoints);
   static ivect make_ghost_zone_vect (CCTK_INT ghost_size,
                       CCTK_INT ghost_size_x, CCTK_INT ghost_size_y,
@@ -413,7 +416,9 @@ namespace Carpet {
       allocate_map (cgh, m, domain_from_coordbase, domain_from_multipatch,
                convergence_factor, refinement_factor,
                base_outerbounds, base_extents,
-               max_refinement_levels, buffer_width, prolongation_order_space,
+               max_refinement_levels, prolongation_order_space,
+               buffer_width,
+               use_outer_buffer_zones, num_integrator_substeps,
                lghosts, ughosts, npoints);
     } 
     
@@ -562,8 +567,12 @@ namespace Carpet {
                CCTK_INT convergence_factor, CCTK_INT refinement_factor,
                CCTK_STRING base_outerbounds, CCTK_STRING base_extents,
                CCTK_INT max_refinement_levels,
-               CCTK_INT buffer_width, CCTK_INT prolongation_order_space,
-               ivect & lghosts, ivect & ughosts, ivect & a_npoints) {
+               CCTK_INT prolongation_order_space,
+               CCTK_INT buffer_width,
+               CCTK_INT use_outer_buffer_zones,
+               CCTK_INT num_integrator_substeps,
+               ivect & lghosts, ivect & ughosts, ivect & a_npoints)
+  {
     // Get boundary description
     jjvect nboundaryzones, is_internal, is_staggered, shiftout;
     int ierr = GetBoundarySpecification (2*dim, &nboundaryzones[0][0],
@@ -795,8 +804,18 @@ namespace Carpet {
                         convergence_factor, vertex_centered, baseext);
     
     // Allocate data hierarchy
+    int const inner_buffer_width = use_outer_buffer_zones ? 0 : buffer_width;
+    ivect const lbuffers
+      = (use_outer_buffer_zones
+         ? lghosts * num_integrator_substeps + buffer_width
+         : 0);
+    ivect const ubuffers
+      = (use_outer_buffer_zones
+         ? ughosts * num_integrator_substeps + buffer_width
+         : 0);
     vdd.at(m) = new dh (*vhh.at(m), lghosts, ughosts,
-                        prolongation_order_space, buffer_width);
+                        prolongation_order_space,
+                        inner_buffer_width, lbuffers, ubuffers);
     
     // Allocate time hierarchy
     vtt.at(m) = new th (*vhh.at(m), timereffacts, 1.0);
@@ -1191,7 +1210,7 @@ namespace Carpet {
                   abaseext);
       
       arrdata.at (group).at (0).dd
-        = new dh (*arrdata.at (group).at (0).hh, alghosts, aughosts, 0, 0);
+        = new dh (*arrdata.at (group).at (0).hh, alghosts, aughosts, 0, 0, 0, 0);
       
       arrdata.at (group).at (0).tt
         = new th (*arrdata.at (group).at (0).hh, atimereffacts, 1.0);
diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc
index 0c5b34e35..16c8f21c5 100644
--- a/Carpet/Carpet/src/modes.cc
+++ b/Carpet/Carpet/src/modes.cc
@@ -95,7 +95,7 @@ namespace Carpet {
 	const ibbox& ext = arrdata.at(group).at(m).dd->boxes.at(ml).at(rl).at(c).exterior;
         
         ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
-          = arrdata.at(group).at(m).dd->lghosts;
+          = arrdata.at(group).at(m).dd->ghosts[0];
         ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
           = base.shape() / base.stride();
         ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
@@ -180,7 +180,7 @@ namespace Carpet {
 //         ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
 //           = deadbeef;
         ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
-          = arrdata.at(group).at(m).dd->lghosts;
+          = arrdata.at(group).at(m).dd->ghosts[0];
         ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
           = deadbeef;
         ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
@@ -298,8 +298,8 @@ namespace Carpet {
     ivect::ref(cgh->cctk_levoff) = (baseext.lower() - coarseext.lower()) / baseext.stride();
     ivect::ref(cgh->cctk_levoffdenom) = 1;
     ivect::ref(cgh->cctk_gsh) = baseext.shape() / baseext.stride();
-    assert (all (vdd.at(map)->lghosts == vdd.at(map)->ughosts));
-    ivect::ref(cgh->cctk_nghostzones) = vdd.at(map)->lghosts;
+    assert (all (vdd.at(map)->ghosts[0] == vdd.at(map)->ghosts[1]));
+    ivect::ref(cgh->cctk_nghostzones) = vdd.at(map)->ghosts[0];
     
     for (int group=0; group<CCTK_NumGroups(); ++group) {
       if (CCTK_GroupTypeI(group) == CCTK_GF) {
@@ -327,7 +327,7 @@ namespace Carpet {
     ivect::ref(cgh->cctk_levoffdenom) = 0;
     ivect::ref(cgh->cctk_gsh) = deadbeef;
 //     ivect::ref(cgh->cctk_nghostzones) = deadbeef;
-    ivect::ref(cgh->cctk_nghostzones) = vdd.at(map)->lghosts;
+    ivect::ref(cgh->cctk_nghostzones) = vdd.at(map)->ghosts[0];
     
     for (int group=0; group<CCTK_NumGroups(); ++group) {
       if (CCTK_GroupTypeI(group) == CCTK_GF) {