Carpet: Implement fully tapered time evolution

Algorithm from Lehner et al.
Original patch concept by Ian Hawke.

darcs-hash:20060603031952-dae7b-897df898e92d7e5807d1df6ab41e3a5d8869deeb.gz

3 files changed, 50 insertions, 11 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index d421b9de8..cb1016fdb 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -389,3 +389,9 @@ BOOLEAN init_3_timelevels "Set up 3 timelevels of initial data" STEERABLE=always
 BOOLEAN adaptive_stepsize "Allow adaptive timestep sizes"
 {
 } "no"
+
+
+
+BOOLEAN use_tapered_grids "Use tapered grids, avoiding time interpolation during evolution"
+{
+} "no"
diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index 7e6489386..121ebf706 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -1,4 +1,6 @@
+#include <algorithm>
 #include <cassert>
+#include <cmath>
 #include <cstdlib>
 
 #include "cctk.h"
@@ -100,7 +102,7 @@ namespace Carpet {
         }
       }
 
-      if (! CCTK_QueryGroupStorageI (cctkGH, group)) {
+      if (not CCTK_QueryGroupStorageI (cctkGH, group)) {
         CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "Cannot synchronise group \"%s\" because it has no storage",
                     groupname);
@@ -115,8 +117,30 @@ namespace Carpet {
 
     if (goodgroups.size() > 0) {
       // prolongate boundaries
-      if (do_prolongate && reflevel > 0) {
-        ProlongateGroupBoundaries (cctkGH, cctk_initial_time, goodgroups);
+      if (reflevel > 0) {
+        if (do_prolongate) {
+          bool local_do_prolongate;
+          if (use_tapered_grids) {
+            // TODO: Check iteration number instead
+            CCTK_REAL mytime;
+            CCTK_REAL parenttime;
+            if (map == -1) {
+              mytime = vtt.at(0)->time (0, reflevel, mglevel);
+              parenttime = vtt.at(0)->time (0, reflevel - 1, mglevel);
+            } else {
+              mytime = vtt.at(map)->time (0, reflevel, mglevel);
+              parenttime = vtt.at(map)->time (0, reflevel - 1, mglevel);
+            }
+            bool const in_sync =
+              abs (mytime - parenttime) < 1.0e-10 * abs (delta_time);
+            local_do_prolongate = in_sync;
+          } else {
+            local_do_prolongate = true;
+          }
+          if (local_do_prolongate) {
+            ProlongateGroupBoundaries (cctkGH, cctk_initial_time, goodgroups);
+          }
+        }
       }
 
       // synchronise ghostzones
@@ -138,7 +162,7 @@ namespace Carpet {
     const CCTK_REAL time
       = (cctkGH->cctk_time - initial_time) / delta_time;
 
-    for (comm_state state; ! state.done(); state.step()) {
+    for (comm_state state; not state.done(); state.step()) {
       for (int group = 0; group < groups.size(); ++group) {
         const int g = groups[group];
         const int grouptype = CCTK_GroupTypeI (g);
@@ -168,7 +192,7 @@ namespace Carpet {
 
     assert (groups.size() > 0);
 
-    for (comm_state state; ! state.done(); state.step()) {
+    for (comm_state state; not state.done(); state.step()) {
       for (int group = 0; group < groups.size(); ++group) {
         const int g = groups[group];
         const int grouptype = CCTK_GroupTypeI (g);
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 520ddb850..f5475f4ff 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -323,6 +323,7 @@ namespace Carpet {
                       CCTK_INT buffer_width,
                       CCTK_INT use_outer_buffer_zones,
                       CCTK_INT num_integrator_substeps,
+                      CCTK_INT use_tapered_grids,
                       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,
@@ -444,6 +445,7 @@ namespace Carpet {
                max_refinement_levels, prolongation_order_space,
                buffer_width,
                use_outer_buffer_zones, num_integrator_substeps,
+               use_tapered_grids,
                lghosts, ughosts, npoints);
     } 
     
@@ -601,6 +603,7 @@ namespace Carpet {
                CCTK_INT buffer_width,
                CCTK_INT use_outer_buffer_zones,
                CCTK_INT num_integrator_substeps,
+               CCTK_INT use_tapered_grids,
                ivect & lghosts, ivect & ughosts, ivect & a_npoints)
   {
     // Get boundary description
@@ -836,13 +839,19 @@ namespace Carpet {
     // 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 - 1) + buffer_width
-         : 0);
+      = (((use_tapered_grids ? refinement_factor : 1)
+          *
+          (use_outer_buffer_zones
+           ? lghosts * num_integrator_substeps + buffer_width
+           : lghosts))
+         - lghosts);
     ivect const ubuffers
-      = (use_outer_buffer_zones
-         ? ughosts * (num_integrator_substeps - 1) + buffer_width
-         : 0);
+      = (((use_tapered_grids ? refinement_factor : 1)
+          *
+          (use_outer_buffer_zones
+           ? ughosts * num_integrator_substeps + buffer_width
+           : ughosts))
+         - ughosts);
     vdd.at(m) = new dh (*vhh.at(m), lghosts, ughosts,
                         prolongation_order_space,
                         inner_buffer_width, lbuffers, ubuffers);