Introduce a parameter to initialise each timelevel separately with

Introduce a parameter to initialise each timelevel separately with
exact initial data.

Do not complain if a grid array has too few timelevels for time
interpolation -- grid arrays are never interpolated in time anyway.

darcs-hash:20030521123023-07bb3-a2e16609a255c931133be9a517363cb1cded29bd.gz

4 files changed, 66 insertions, 22 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index b8f7b5a79..8a9b563f2 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn Carpet
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.28 2003/05/02 14:22:31 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.29 2003/05/21 14:30:23 schnetter Exp $
 
 shares: Cactus
 
@@ -229,6 +229,10 @@ BOOLEAN periodic_z "do not use this parameter"
 
 private:
 
+BOOLEAN init_each_timelevel "Call initial data routines once for each timelevel"
+{
+} "no"
+
 BOOLEAN init_3_timelevels "Set up 3 timelevels of initial data"
 {
 } "no"
diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index e6fb610a2..c3f59803b 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -10,7 +10,7 @@
 #include "carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.15 2003/01/03 14:11:56 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.16 2003/05/21 14:30:24 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_Carpet_Comm_cc);
 }
 
@@ -62,7 +62,7 @@ namespace Carpet {
               if (false) {
                 const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel);
                 const CCTK_REAL time2 = cgh->cctk_time / base_delta_time;
-                assert (fabs(time1 - time2) < 1e-10);
+                assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(base_delta_time))) < 1e-12);
               }
               
               arrdata[group].data[var]->ref_bnd_prolongate
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 266ab8958..a8c1f2f73 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -12,7 +12,7 @@
 #include "carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.27 2003/05/13 16:31:47 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.28 2003/05/21 14:30:24 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_Carpet_Initialise_cc);
 }
 
@@ -92,11 +92,48 @@ namespace Carpet {
 	  cgh->cctk_delta_time = base_delta_time / reflevelfact * mglevelfact;
 	}
 	
-	// Set up the initial data
-	Waypoint ("%*sScheduling INITIAL", 2*reflevel, "");
-	CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction);
-	Waypoint ("%*sScheduling POSTINITIAL", 2*reflevel, "");
-	CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction);
+        if (! init_each_timelevel) {
+          
+          // Set up the initial data
+          Waypoint ("%*sScheduling INITIAL", 2*reflevel, "");
+          CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction);
+          Waypoint ("%*sScheduling POSTINITIAL", 2*reflevel, "");
+          CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction);
+          
+        } else {
+          // init_each_timelevel
+          
+	  tt->set_delta
+            (reflevel, mglevel, - tt->get_delta (reflevel, mglevel));
+	  tt->advance_time (reflevel, mglevel);
+	  tt->advance_time (reflevel, mglevel);
+	  tt->advance_time (reflevel, mglevel);
+	  tt->set_delta
+            (reflevel, mglevel, - tt->get_delta (reflevel, mglevel));
+          
+          for (int tl=-2; tl<=0; ++tl) {
+            
+            // Advance level times
+            tt->advance_time (reflevel, mglevel);
+            cgh->cctk_time
+              = tt->time (0, reflevel, mglevel) * base_delta_time;
+            
+            // Cycle time levels
+            CycleTimeLevels (cgh);
+            
+            Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+                      cgh->cctk_time,
+                      do_global_mode ? "   (global time)" : "");
+            
+            // Set up the initial data
+            Waypoint ("%*sScheduling INITIAL", 2*reflevel, "");
+            CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction);
+            Waypoint ("%*sScheduling POSTINITIAL", 2*reflevel, "");
+            CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction);
+            
+          } // for tl
+          
+        } // init_each_timelevel
 	
 	// Poststep
 	Waypoint ("%*sScheduling POSTSTEP", 2*reflevel, "");
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index 1b4ed7ecf..c17842da5 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -10,7 +10,7 @@
 #include "carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.17 2003/05/15 16:31:46 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.18 2003/05/21 14:30:25 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_Carpet_Storage_cc);
 }
 
@@ -39,18 +39,19 @@ namespace Carpet {
     const int group = CCTK_GroupIndex(groupname);
     assert (group>=0 && group<CCTK_NumGroups());
     
+    const int grouptype = CCTK_GroupTypeI(group);
+    
     // No storage change in local mode
-    assert (! (component!=-1 && CCTK_GroupTypeI(group)==CCTK_GF));
+    assert (! (component!=-1 && grouptype==CCTK_GF));
     
     if (CCTK_QueryGroupStorageI(cgh, group)) {
       // storage was enabled previously
       return 1;
     }
-    
+
     // Check whether this group has transfer operators
-    if (! arrdata[group].do_transfer) {
-      const int grouptype = CCTK_GroupTypeI(group);
-      if (grouptype == CCTK_GF) {
+    if (grouptype == CCTK_GF) {
+      if (! arrdata[group].do_transfer) {
         const int var = CCTK_FirstVarIndexI(group);
         const int vartype = CCTK_VarTypeI(var);
         const char * vartypename = CCTK_VarTypeName(vartype);
@@ -74,12 +75,14 @@ namespace Carpet {
     const int my_prolongation_order_time
       = num_tl==1 ? 0 : prolongation_order_time;
     
-    if (max_refinement_levels > 1) {
-      if (num_tl <= my_prolongation_order_time) {
-	CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "There are not enough time levels for the desired temporal prolongation order in the grid function group \"%s\".  With Carpet::prolongation_order_time=%d, you need at least %d time levels.",
-		    CCTK_GroupName(group),
-		    prolongation_order_time, prolongation_order_time+1);
+    if (grouptype == CCTK_GF) {
+      if (max_refinement_levels > 1) {
+        if (num_tl <= my_prolongation_order_time) {
+          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "There are not enough time levels for the desired temporal prolongation order in the grid function group \"%s\".  With Carpet::prolongation_order_time=%d, you need at least %d time levels.",
+                      CCTK_GroupName(group),
+                      prolongation_order_time, prolongation_order_time+1);
+        }
       }
     }
     
@@ -100,7 +103,7 @@ namespace Carpet {
 	UnsupportedVarType(n);
       } // switch
       
-      if (CCTK_GroupTypeI(group) != CCTK_GF) {
+      if (grouptype != CCTK_GF) {
         for (int tl=0; tl<num_tl; ++tl) {
           int const c = CCTK_MyProc(cgh);
           cgh->data[n][tl] = ((*arrdata[group].data[var]) (-tl, 0, c, 0)->storage());