Carpet: Correct problems with negative time steps

When converting from physicsl time steps (delta_time) to Carpet's time
steps (th::delta_time), do not divide by abs(delta_time), but by
delta_time instead.

When evolving backwards in time while initialising three time levels,
change the sign of the physical time steps only, not the sign of
Carpet's time steps.

darcs-hash:20050809084940-891bb-2517b1568696c71278a98db6261515817a90247a.gz

3 files changed, 13 insertions, 19 deletions

diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index 12d7e1394..0e0965ae4 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -1,5 +1,4 @@
 #include <cassert>
-#include <cmath>
 #include <cstdlib>
 
 #include "cctk.h"
@@ -140,7 +139,7 @@ namespace Carpet {
 
     // use the current time here (which may be modified by the user)
     const CCTK_REAL time
-      = (cctkGH->cctk_time - initial_time) / abs(delta_time);
+      = (cctkGH->cctk_time - initial_time) / delta_time;
 
     for (comm_state state(groups.vartype); ! state.done(); state.step()) {
       for (int group = 0; group < groups.members.size(); ++group) {
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 4c2adca70..ade713ecc 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -108,8 +108,6 @@ 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);
       }
     }
@@ -448,8 +446,13 @@ namespace Carpet {
   // Use Scott Hawley's algorithm to get two extra timelevels of data
   void initialise_3tl (cGH * const cgh)
   {
+    DECLARE_CCTK_PARAMETERS;
+    
     Waypoint ("Initialising three timelevels");
 
+    // TODO: ensure that there are 3 timelevels
+    assert (prolongation_order_time == 2);
+
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
@@ -502,10 +505,11 @@ namespace Carpet {
   {
     Waypoint ("Advancing time");
     
+    cgh->cctk_time
+      = global_time + delta_time * mglevelfact / timereflevelfact;
     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);
   }
@@ -539,12 +543,6 @@ namespace Carpet {
     BEGIN_MGLEVEL_LOOP(cgh) {
       BEGIN_REFLEVEL_LOOP (cgh) {
 
-        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 + delta_time * mglevelfact / timereflevelfact;
 
@@ -603,11 +601,11 @@ namespace Carpet {
   {
     Waypoint ("Advancing time");
     
+    cgh->cctk_time
+      = global_time + 2 * delta_time * mglevelfact / timereflevelfact;
     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);
   }
@@ -644,10 +642,7 @@ namespace Carpet {
       BEGIN_REFLEVEL_LOOP (cgh) {
 
         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_time (reflevel, mglevel, 0);
         }
         cgh->cctk_time = global_time;
 
diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc
index 78a3857b1..42d4d8720 100644
--- a/Carpet/Carpet/src/Restrict.cc
+++ b/Carpet/Carpet/src/Restrict.cc
@@ -92,8 +92,8 @@ namespace Carpet {
 
           const CCTK_REAL time1 = vtt.at(m)->time (0, reflevel, mglevel);
           const CCTK_REAL time2
-            = (cgh->cctk_time - cctk_initial_time) / fabs(delta_time);
-          assert (fabs(time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time)) < 1e-12);
+            = (cgh->cctk_time - cctk_initial_time) / delta_time;
+          assert (abs(time1 - time2) / (abs(time1) + abs(time2) + abs(cgh->cctk_delta_time)) < 1e-12);
 
           for (int v = 0; v < arrdata.at(group).at(m).data.size(); ++v) {
             ggf *const gv = arrdata.at(group).at(m).data.at(v);