Carpet: Make 3-timelevel-initialisation work again

3-timelevel-initialisation was broken: the refinement levels were
initialised in the wrong order, so that the method was formally only
first order convergent.

Reorder how the individual steps of the algorithm traverse the
refinement levels.

darcs-hash:20050725212427-891bb-72d0e8554e98e7126187b61ef87fe447085dc4af.gz

2 files changed, 149 insertions, 132 deletions

diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 70cdf0e58..4c2adca70 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -34,10 +34,13 @@ namespace Carpet {
   static void initialise_III (cGH * cgh, int rl);
 
   static void initialise_3tl (cGH * cgh);
+  static void initialise_3tl_advance_time (cGH * cgh);
   static void initialise_3tl_evolve_Ia (cGH * cgh);
   static void initialise_3tl_flip_timelevels (cGH * cgh);
   static void initialise_3tl_evolve_Ib (cGH * cgh);
-  static void initialise_3tl_evolve_IIb_Ic (cGH * cgh);
+  static void initialise_3tl_evolve_IIb (cGH * cgh);
+  static void initialise_3tl_advance_time_2 (cGH * cgh);
+  static void initialise_3tl_evolve_Ic (cGH * cgh);
   static void initialise_3tl_flip_timelevels_back (cGH * cgh);
 
   static void print_internal_data ();
@@ -105,6 +108,8 @@ namespace Carpet {
       print_internal_data ();
 
       if (init_3_timelevels) {
+#warning "TODO: ensure that there are 3 timelevels"
+#warning "TODO: ensure that prolongation_order_time = 2"
         initialise_3tl (cgh);
       }
     }
@@ -445,199 +450,211 @@ namespace Carpet {
   {
     Waypoint ("Initialising three timelevels");
 
-    initialise_3tl_evolve_Ia (cgh);
-
-    delta_time *= -1;
-    initialise_3tl_flip_timelevels (cgh);
-
-    initialise_3tl_evolve_Ib (cgh);
-
-    Waypoint ("Hourglass structure in place");
-
-    initialise_3tl_evolve_IIb_Ic (cgh);
-
-    delta_time *= -1;
-    initialise_3tl_flip_timelevels_back (cgh);
-
-    Waypoint ("Finished initialising three timelevels");
-  }
-
-  void initialise_3tl_evolve_Ia (cGH * const cgh)
-  {
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==0;
         do_meta_mode = do_global_mode and mglevel==mglevels-1;
 
-        // Advance times
-        for (int m=0; m<maps; ++m) {
-          vtt.at(m)->advance_time (reflevel, mglevel);
-        }
-        cgh->cctk_time
-          = global_time + delta_time * mglevelfact / timereflevelfact;
-        CycleTimeLevels (cgh);
-
-        Waypoint ("Initialisation 3TL evolution I (a) (forwards) at iteration"
-                  " %d time %g%s%s",
-                  cgh->cctk_iteration, (double)cgh->cctk_time,
-                  (do_global_mode ? " (global)" : ""),
-                  (do_meta_mode ? " (meta)" : ""));
-
-        CalculateChecksums (cgh, allbutcurrenttime);
-        Poison (cgh, currenttimebutnotifonly);
-
-        // Evolve forward
-        Checkpoint ("Scheduling PRESTEP");
-        CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
-        Checkpoint ("Scheduling EVOL");
-        CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
-
-        PoisonCheck (cgh, currenttime);
+        initialise_3tl_advance_time (cgh);
+        initialise_3tl_evolve_Ia (cgh);
 
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
-  }
 
-  void initialise_3tl_flip_timelevels (cGH * const cgh)
-  {
+    initialise_3tl_flip_timelevels (cgh);
+
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==0;
         do_meta_mode = do_global_mode and mglevel==mglevels-1;
 
-        // Flip time levels
-        Waypoint ("Flipping timelevels");
-        FlipTimeLevels (cgh);
-
-        cgh->cctk_time
-          = global_time + delta_time * mglevelfact / timereflevelfact;
+        initialise_3tl_evolve_Ib (cgh);
 
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
-  }
 
-  void initialise_3tl_evolve_Ib (cGH * const cgh)
-  {
-    for (int rl=0; rl<reflevels; ++rl) {
+    Waypoint ("Hourglass structure in place");
+
+    for (int rl=reflevels-1; rl>=0; --rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==0;
         do_meta_mode = do_global_mode and mglevel==mglevels-1;
 
-        Waypoint ("Initialisation 3TL evolution I (b) (backwards) at iteration"
-                  " %d time %g%s%s",
-                  cgh->cctk_iteration, (double)cgh->cctk_time,
-                  (do_global_mode ? " (global)" : ""),
-                  (do_meta_mode ? " (meta)" : ""));
-
-        // Checking
-        CalculateChecksums (cgh, allbutcurrenttime);
-        Poison (cgh, currenttimebutnotifonly);
-
-        // Evolve backward
-        Checkpoint ("Scheduling PRESTEP");
-        CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
-        Checkpoint ("Scheduling EVOL");
-        CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
-
-        // Checking
-        PoisonCheck (cgh, alltimes);
+        initialise_3tl_evolve_IIb (cgh);
+        initialise_3tl_advance_time_2 (cgh);
+        initialise_3tl_evolve_Ic (cgh);
 
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
+
+    initialise_3tl_flip_timelevels_back (cgh);
+
+    Waypoint ("Finished initialising three timelevels");
   }
 
-  // Evolve each level "backwards" one more timestep
-  // Starting with the finest level and proceeding to the coarsest
-  void initialise_3tl_evolve_IIb_Ic (cGH * const cgh)
+  void initialise_3tl_advance_time (cGH * const cgh)
   {
-    for (int rl=reflevels-1; rl>=0; --rl) {
-      BEGIN_MGLEVEL_LOOP(cgh) {
-        enter_level_mode (cgh, rl);
-        do_global_mode = reflevel==0;
-        do_meta_mode = do_global_mode and mglevel==mglevels-1;
-
-        Waypoint ("Initialisation 3TL evolution II (b) (backwards) at iteration"
-                  " %d time %g%s%s",
-                  cgh->cctk_iteration, (double)cgh->cctk_time,
-                  (do_global_mode ? " (global)" : ""),
-                  (do_meta_mode ? " (meta)" : ""));
+    Waypoint ("Advancing time");
+    
+    for (int m=0; m<maps; ++m) {
+      vtt.at(m)->advance_time (reflevel, mglevel);
+    }
+    cgh->cctk_time = global_time + delta_time * mglevelfact / timereflevelfact;
+    
+    CycleTimeLevels (cgh);
+  }
 
-        Restrict (cgh);
+  void initialise_3tl_evolve_Ia (cGH * const cgh)
+  {
+    Waypoint ("Initialisation 3TL evolution I (a) (forwards) at iteration"
+              " %d time %g%s%s",
+              cgh->cctk_iteration, (double)cgh->cctk_time,
+              (do_global_mode ? " (global)" : ""),
+              (do_meta_mode ? " (meta)" : ""));
+
+    CalculateChecksums (cgh, allbutcurrenttime);
+    Poison (cgh, currenttimebutnotifonly);
+
+    // Evolve forward
+    Checkpoint ("Scheduling PRESTEP");
+    CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
+    Checkpoint ("Scheduling EVOL");
+    CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
+
+    PoisonCheck (cgh, currenttime);
+  }
 
-        if (rl < reflevels-1) {
-          Checkpoint ("Scheduling POSTRESTRICT");
-          CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction);
-        }
+  void initialise_3tl_flip_timelevels (cGH * const cgh)
+  {
+    Waypoint ("Flipping timelevels");
 
-        Checkpoint ("Scheduling POSTSTEP");
-        CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
+    delta_time *= -1;
 
-        PoisonCheck (cgh, alltimes);
+    BEGIN_MGLEVEL_LOOP(cgh) {
+      BEGIN_REFLEVEL_LOOP (cgh) {
 
-        // Advance times
         for (int m=0; m<maps; ++m) {
+          vtt.at(m)->set_delta
+            (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
+          vtt.at(m)->advance_time (reflevel, mglevel);
           vtt.at(m)->advance_time (reflevel, mglevel);
         }
         cgh->cctk_time
-          = global_time + 2 * delta_time * mglevelfact / timereflevelfact;
-        CycleTimeLevels (cgh);
+          = global_time + delta_time * mglevelfact / timereflevelfact;
 
-        Waypoint ("Initialisation 3TL evolution I (c) (backwards) at iteration"
-                  " %d time %g%s%s",
-                  cgh->cctk_iteration, (double)cgh->cctk_time,
-                  (do_global_mode ? " (global)" : ""),
-                  (do_meta_mode ? " (meta)" : ""));
+        FlipTimeLevels (cgh);
 
-        CalculateChecksums (cgh, allbutcurrenttime);
-        Poison (cgh, currenttimebutnotifonly);
+      } END_REFLEVEL_LOOP;
+    } END_MGLEVEL_LOOP;
+  }
 
-        // Evolve backward
-        Checkpoint ("Scheduling PRESTEP");
-        CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
-        Checkpoint ("Scheduling EVOL");
-        CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
-        Checkpoint ("Scheduling POSTSTEP");
-        CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
+  void initialise_3tl_evolve_Ib (cGH * const cgh)
+  {
+    Waypoint ("Initialisation 3TL evolution I (b) (backwards) at iteration"
+              " %d time %g%s%s",
+              cgh->cctk_iteration, (double)cgh->cctk_time,
+              (do_global_mode ? " (global)" : ""),
+              (do_meta_mode ? " (meta)" : ""));
+
+    // Checking
+    CalculateChecksums (cgh, allbutcurrenttime);
+    Poison (cgh, currenttimebutnotifonly);
+
+    // Evolve backward
+    Checkpoint ("Scheduling PRESTEP");
+    CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
+    Checkpoint ("Scheduling EVOL");
+    CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
+
+    // Checking
+    PoisonCheck (cgh, alltimes);
+  }
 
-        PoisonCheck (cgh, alltimes);
+  // Evolve backwards one more timestep
+  // Starting with the finest level and proceeding to the coarsest
+  void initialise_3tl_evolve_IIb (cGH * const cgh)
+  {
+    Waypoint ("Initialisation 3TL evolution II (b) (backwards) at iteration"
+              " %d time %g%s%s",
+              cgh->cctk_iteration, (double)cgh->cctk_time,
+              (do_global_mode ? " (global)" : ""),
+              (do_meta_mode ? " (meta)" : ""));
 
-        leave_level_mode (cgh);
-      } END_MGLEVEL_LOOP;
+    Restrict (cgh);
+
+    if (reflevel < reflevels-1) {
+      Checkpoint ("Scheduling POSTRESTRICT");
+      CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction);
+    }
+
+    Checkpoint ("Scheduling POSTSTEP");
+    CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
+
+    PoisonCheck (cgh, alltimes);
+  }
+
+  void initialise_3tl_advance_time_2 (cGH * const cgh)
+  {
+    Waypoint ("Advancing time");
+    
+    for (int m=0; m<maps; ++m) {
+      vtt.at(m)->advance_time (reflevel, mglevel);
     }
+    cgh->cctk_time
+      = global_time + 2 * delta_time * mglevelfact / timereflevelfact;
+    
+    CycleTimeLevels (cgh);
   }
 
+  void initialise_3tl_evolve_Ic (cGH * const cgh)
+  {
+    Waypoint ("Initialisation 3TL evolution I (c) (backwards) at iteration"
+              " %d time %g%s%s",
+              cgh->cctk_iteration, (double)cgh->cctk_time,
+              (do_global_mode ? " (global)" : ""),
+              (do_meta_mode ? " (meta)" : ""));
+
+    CalculateChecksums (cgh, allbutcurrenttime);
+    Poison (cgh, currenttimebutnotifonly);
+
+    // Evolve backward
+    Checkpoint ("Scheduling PRESTEP");
+    CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
+    Checkpoint ("Scheduling EVOL");
+    CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
+    Checkpoint ("Scheduling POSTSTEP");
+    CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
+
+    PoisonCheck (cgh, alltimes);
+  }
+  
   void initialise_3tl_flip_timelevels_back (cGH * const cgh)
   {
-    for (int rl=0; rl<reflevels; ++rl) {
-      BEGIN_MGLEVEL_LOOP(cgh) {
-        enter_level_mode (cgh, rl);
-        do_global_mode = reflevel==0;
-        do_meta_mode = do_global_mode and mglevel==mglevels-1;
+    Waypoint ("Flipping timelevels back");
 
-        // Flip time levels back
-        Waypoint ("Flipping timelevels back");
-        FlipTimeLevels (cgh);
+    delta_time *= -1;
+
+    BEGIN_MGLEVEL_LOOP(cgh) {
+      BEGIN_REFLEVEL_LOOP (cgh) {
 
-        // Invert level times back
         for (int m=0; m<maps; ++m) {
           vtt.at(m)->set_delta
             (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
           vtt.at(m)->advance_time (reflevel, mglevel);
           vtt.at(m)->advance_time (reflevel, mglevel);
-          vtt.at(m)->set_delta
-            (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
         }
         cgh->cctk_time = global_time;
 
-        leave_level_mode (cgh);
-      } END_MGLEVEL_LOOP;
-    }
+        FlipTimeLevels (cgh);
+
+      } END_REFLEVEL_LOOP;
+    } END_MGLEVEL_LOOP;
   }
 
 
diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc
index d4d1e962e..78a3857b1 100644
--- a/Carpet/Carpet/src/Restrict.cc
+++ b/Carpet/Carpet/src/Restrict.cc
@@ -92,7 +92,7 @@ namespace Carpet {
 
           const CCTK_REAL time1 = vtt.at(m)->time (0, reflevel, mglevel);
           const CCTK_REAL time2
-            = (cgh->cctk_time - cctk_initial_time) / delta_time;
+            = (cgh->cctk_time - cctk_initial_time) / fabs(delta_time);
           assert (fabs(time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time)) < 1e-12);
 
           for (int v = 0; v < arrdata.at(group).at(m).data.size(); ++v) {