call initial data checkpoint, analysis and output routines only for refinement levels which are to be updated at the current iteration

This patch effects only initial data after recovery: it prevents the analysis
and output of refinement levels which weren't updated at the iteration when
the checkpoint was created.

darcs-hash:20041209135038-3fd61-96f33595a27f7705734ae1822427bbf02e7fe389.gz

1 files changed, 145 insertions, 136 deletions

diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index d2c7e19e0..8d6003ee9 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -19,16 +19,16 @@ extern "C" {
 
 
 namespace Carpet {
-  
+
   using namespace std;
-  
-  
+
+
   static void output_the_grid_structure ( cGH* cgh );
   static void register_coordinates_and_check_parameters ( cGH* cgh );
   static void recovery_I ( cGH* cgh, int rl );
   static void recovery_Regrid ( cGH* cgh, int rl );
   static void recovery_II ( cGH* cgh );
-  static void 
+  static void
   initialisation_I ( cGH* cgh, int rl, int init_each_timelevel );
   static void initialise_rewind ( cGH* cgh, int num_tl );
   static void initialise_Schedule_INITIAL ( cGH* cgh, int num_tl );
@@ -43,64 +43,64 @@ namespace Carpet {
   static void initialise_evolve_3TL_backwards_IIb_Ic ( cGH* cgh );
   static void initialise_Flip_Timelevels_back ( cGH* cgh );
   static void initialisation_III ( cGH* cgh );
-  
+
   int Initialise (tFleshConfig* fc)
   {
     DECLARE_CCTK_PARAMETERS;
-    
+
     const int convlev = 0;
     cGH* const cgh = CCTK_SetupGH (fc, convlev);
     CCTKi_AddGH (fc, convlev, cgh);
-    
+
     do_global_mode = true;
     do_meta_mode = true;
     global_time = cctk_initial_time;
     delta_time = 1.0;
-    
+
     cgh->cctk_iteration = 0;
     cgh->cctk_time = global_time;
     cgh->cctk_delta_time = delta_time;
 
     // Delay checkpoint until MPI has been initialised
     Waypoint ("Starting initialisation");
-    
+
     CCTKi_ScheduleGHInit (cgh); // Enable storage and communication
-    
+
     CCTKi_InitGHExtensions (cgh);
-    
+
     output_the_grid_structure (cgh);
-    
+
     register_coordinates_and_check_parameters (cgh);
-    
+
     if (fc->recovered) {
-      
+
       for (int rl=0; rl<reflevels; ++rl) {
         recovery_I (cgh, rl);
         recovery_Regrid (cgh, rl);
-      } 
+      }
 
       recovery_II (cgh);
-      
+
     } else {
 
       for (int rl=0; rl<reflevels; ++rl) {
         initialisation_I (cgh, rl, init_each_timelevel);
         initialise_Regrid (cgh, rl, prolongate_initial_data);
       }
-      
+
       initialise_Restrict (cgh);
-      
+
       initialisation_II (cgh);
-    
+
       if (init_3_timelevels) {
         get_two_extra_timelevels_of_data (cgh);
-      } 
+      }
     }
-    
+
     initialisation_III (cgh);
- 
+
     Waypoint ("Done with initialisation");
-    
+
     return 0;
   }
 
@@ -121,143 +121,143 @@ namespace Carpet {
     BEGIN_MGLEVEL_LOOP(cgh) {
       do_global_mode = true;
       do_meta_mode = mglevel==mglevels-1;
-      
+
       // Register coordinates
       Checkpoint ("Scheduling CCTK_WRAGH");
       CCTK_ScheduleTraverse ("CCTK_WRAGH", cgh, CallFunction);
-      
+
       // Check parameters
       Checkpoint ("Scheduling PARAMCHECK");
       CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction);
       CCTKi_FinaliseParamWarn();
     } END_MGLEVEL_LOOP;
   }
-    
+
   void recovery_I (cGH* cgh, int rl)
-  { 
+  {
     BEGIN_MGLEVEL_LOOP(cgh) {
       enter_level_mode (cgh, rl);
       do_global_mode = reflevel==0;
       do_meta_mode = do_global_mode && mglevel==mglevels-1;
-      
+
       cgh->cctk_time = global_time;
-      
+
       Waypoint ("Recovering I at iteration %d time %g%s%s",
                 cgh->cctk_iteration, (double)cgh->cctk_time,
                 (do_global_mode ? " (global)" : ""),
                 (do_meta_mode ? " (meta)" : ""));
-      
+
       // Set up the grids
       Checkpoint ("Scheduling BASEGRID");
       CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction);
-          
+
       // Recover
       Checkpoint ("Scheduling RECOVER_VARIABLES");
       CCTK_ScheduleTraverse ("CCTK_RECOVER_VARIABLES", cgh, CallFunction);
-      
+
       leave_level_mode (cgh);
     } END_MGLEVEL_LOOP;
   }
 
   void recovery_Regrid (cGH* cgh, int rl)
-  { 
+  {
     bool did_regrid = false;
     {
       const int ml=0;
       enter_global_mode (cgh, ml);
       enter_level_mode (cgh, rl);
-      
+
       // Regrid
       Checkpoint ("Regrid");
       did_regrid |= Regrid (cgh, true, false);
-      
+
       leave_level_mode (cgh);
       leave_global_mode (cgh);
     } // ml
-    
+
     if (did_regrid) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = true;
         do_meta_mode = do_global_mode && mglevel==mglevels-1;
-        
+
         Waypoint ("Postregrid at iteration %d time %g%s%s",
                   cgh->cctk_iteration, (double)cgh->cctk_time,
                   (do_global_mode ? " (global)" : ""),
                   (do_meta_mode ? " (meta)" : ""));
-        
+
         // Postregrid
         Checkpoint ("Scheduling POSTREGRID");
         CCTK_ScheduleTraverse ("CCTK_POSTREGRID", cgh, CallFunction);
-        
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     } // if did_regrid
   }
 
   void recovery_II (cGH* cgh)
-  { 
+  {
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==reflevels-1;
         do_meta_mode = do_global_mode && mglevel==mglevels-1;
-      
+
         Waypoint ("Recovering II at iteration %d time %g%s%s",
                   cgh->cctk_iteration, (double)cgh->cctk_time,
                   (do_global_mode ? " (global)" : ""),
                   (do_meta_mode ? " (meta)" : ""));
-      
+
         // Post recover
         Checkpoint ("Scheduling POST_RECOVER_VARIABLES");
         CCTK_ScheduleTraverse
           ("CCTK_POST_RECOVER_VARIABLES", cgh, CallFunction);
-      
+
         // Checking
         PoisonCheck (cgh, alltimes);
         CheckChecksums (cgh, allbutcurrenttime);
-      
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
-  
+
   void initialisation_I (cGH* cgh, int rl, int init_each_timelevel)
-  { 
+  {
     BEGIN_MGLEVEL_LOOP(cgh) {
       enter_level_mode (cgh, rl);
       do_global_mode = reflevel==0;
       do_meta_mode = do_global_mode && mglevel==mglevels-1;
-      
+
       cgh->cctk_time = global_time;
-      
+
       Waypoint ("Initialisation I at iteration %d time %g%s%s",
                 cgh->cctk_iteration, (double)cgh->cctk_time,
                 (do_global_mode ? " (global)" : ""),
                 (do_meta_mode ? " (meta)" : ""));
-      
+
       // Checking
       Poison (cgh, alltimes);
-      
+
       // Set up the grids
       Checkpoint ("Scheduling BASEGRID");
       CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction);
-      
+
       const int num_tl = init_each_timelevel ? 3 : 1;
-      
+
       initialise_rewind (cgh, num_tl);
 
       initialise_Schedule_INITIAL (cgh, num_tl);
-      
+
       // Checking
       PoisonCheck (cgh, currenttime);
-      
+
       leave_level_mode (cgh);
     } END_MGLEVEL_LOOP;
   }
-  
+
   void initialise_rewind (cGH* cgh, int num_tl)
-  { 
+  {
     for (int m=0; m<maps; ++m) {
       vtt.at(m)->set_delta
         (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
@@ -270,14 +270,14 @@ namespace Carpet {
         (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
       FlipTimeLevels (cgh);
     }
-  } 
-  
+  }
+
   void initialise_Schedule_INITIAL (cGH* cgh, int num_tl)
-  { 
+  {
     const bool outer_do_global_mode = do_global_mode;
     for (int tl=num_tl-1; tl>=0; --tl) {
       do_global_mode = outer_do_global_mode && tl==0;
-      
+
       // Advance times
       for (int m=0; m<maps; ++m) {
         vtt.at(m)->advance_time (reflevel, mglevel);
@@ -285,93 +285,93 @@ namespace Carpet {
       cgh->cctk_time
         = global_time - tl * delta_time * mglevelfact / reflevelfact;
       CycleTimeLevels (cgh);
-      
+
       // Set up the initial data
       Checkpoint ("Scheduling INITIAL");
       CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction);
-      
+
     } // for tl
     do_global_mode = outer_do_global_mode;
-  } 
+  }
 
   void initialise_Regrid (cGH* cgh, int rl, int prolongate_initial_data)
-  { 
+  {
     bool did_regrid = false;
     {
       const int ml=0;
       enter_global_mode (cgh, ml);
       enter_level_mode (cgh, rl);
-      
+
       // Regrid
       Checkpoint ("Regrid");
       did_regrid |= Regrid (cgh, false, prolongate_initial_data);
-      
+
       leave_level_mode (cgh);
       leave_global_mode (cgh);
     } // ml
-    
+
     if (did_regrid) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = true;
         do_meta_mode = do_global_mode && mglevel==mglevels-1;
-        
+
         Waypoint ("Postregrid at iteration %d time %g%s%s",
                   cgh->cctk_iteration, (double)cgh->cctk_time,
                   (do_global_mode ? " (global)" : ""),
                   (do_meta_mode ? " (meta)" : ""));
-        
+
         // Postregrid
         Checkpoint ("Scheduling POSTREGRID");
         CCTK_ScheduleTraverse ("CCTK_POSTREGRID", cgh, CallFunction);
-        
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     } // if did_regrid
   }
 
   void initialise_Restrict (cGH* cgh)
-  { 
+  {
     for (int rl=reflevels-1; rl>=0; --rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
-        
+
         Waypoint ("Initialisation/Restrict at iteration %d time %g",
                   cgh->cctk_iteration, (double)cgh->cctk_time);
-        
+
         Restrict (cgh);
-        
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
 
   void initialisation_II (cGH* cgh)
-  { 
+  {
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==reflevels-1;
         do_meta_mode = do_global_mode && mglevel==mglevels-1;
-      
+
         Waypoint ("Initialisation II at iteration %d time %g%s%s",
                   cgh->cctk_iteration, (double)cgh->cctk_time,
                   (do_global_mode ? " (global)" : ""),
                   (do_meta_mode ? " (meta)" : ""));
-      
+
         Checkpoint ("Scheduling POSTRESTRICTINITIAL");
         CCTK_ScheduleTraverse
           ("CCTK_POSTRESTRICTINITIAL", cgh, CallFunction);
-      
+
         Checkpoint ("Scheduling POSTINITIAL");
         CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction);
-      
+
         Checkpoint ("Scheduling POSTSTEP");
         CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
-      
+
         PoisonCheck (cgh, alltimes);
         CheckChecksums (cgh, allbutcurrenttime);
-      
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
@@ -380,34 +380,34 @@ namespace Carpet {
   // Use Scott Hawley's algorithm for getting two extra
   // timelevels of data
   void get_two_extra_timelevels_of_data (cGH* cgh)
-  { 
+  {
         Waypoint ("Initialising three timelevels");
-      
+
         initialise_3_Timelevels (cgh);
-      
+
         delta_time *= -1;
         initialise_Flip_Timelevels (cgh);
-      
+
         initialise_evolve_3TL_backwards_Ib (cgh);
-      
+
         Waypoint ("Hourglass structure in place");
-      
+
         initialise_evolve_3TL_backwards_IIb_Ic (cgh);
 
         delta_time *= -1;
         initialise_Flip_Timelevels_back (cgh);
-      
+
         Waypoint ("Finished initialising three timelevels");
   }
 
   void initialise_3_Timelevels (cGH* 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 && mglevel==mglevels-1;
-    
+
         // Advance times
         for (int m=0; m<maps; ++m) {
           vtt.at(m)->advance_time (reflevel, mglevel);
@@ -415,76 +415,76 @@ namespace Carpet {
         cgh->cctk_time
           = global_time + delta_time * mglevelfact / reflevelfact;
         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);
-    
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
 
   void initialise_Flip_Timelevels (cGH* 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 && mglevel==mglevels-1;
-    
+
         // Flip time levels
         Waypoint ("Flipping timelevels");
         FlipTimeLevels (cgh);
-    
+
         cgh->cctk_time
           = global_time + delta_time * mglevelfact / reflevelfact;
-    
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
 
   void initialise_evolve_3TL_backwards_Ib (cGH* 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 && 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);
-    
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
@@ -493,29 +493,29 @@ namespace Carpet {
   // Evolve each level "backwards" one more timestep
   // Starting with the finest level and proceeding to the coarsest
   void initialise_evolve_3TL_backwards_IIb_Ic (cGH* 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 && 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)" : ""));
-      
+
         Restrict (cgh);
-      
+
         Checkpoint ("Scheduling POSTRESTRICT");
         CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction);
-      
+
         Checkpoint ("Scheduling POSTSTEP");
         CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
-    
+
         PoisonCheck (cgh, alltimes);
-    
+
         // Advance times
         for (int m=0; m<maps; ++m) {
           vtt.at(m)->advance_time (reflevel, mglevel);
@@ -523,16 +523,16 @@ namespace Carpet {
         cgh->cctk_time
           = global_time + 2 * delta_time * mglevelfact / reflevelfact;
         CycleTimeLevels (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);
@@ -540,26 +540,26 @@ namespace Carpet {
         CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction);
         Checkpoint ("Scheduling POSTSTEP");
         CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
-    
+
         PoisonCheck (cgh, alltimes);
-    
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
 
   void initialise_Flip_Timelevels_back (cGH* 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 && mglevel==mglevels-1;
-    
+
         // Flip time levels back
         Waypoint ("Flipping timelevels back");
         FlipTimeLevels (cgh);
-    
+
         // Invert level times back
         for (int m=0; m<maps; ++m) {
           vtt.at(m)->set_delta
@@ -570,37 +570,46 @@ namespace Carpet {
             (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
         }
         cgh->cctk_time = global_time;
-    
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }
   }
 
   void initialisation_III (cGH* cgh)
-  { 
+  {
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cgh) {
         enter_level_mode (cgh, rl);
         do_global_mode = reflevel==reflevels-1;
         do_meta_mode = do_global_mode && mglevel==mglevels-1;
-      
+
         Waypoint ("Initialisation III at iteration %d time %g%s%s",
                   cgh->cctk_iteration, (double)cgh->cctk_time,
                   (do_global_mode ? " (global)" : ""),
                   (do_meta_mode ? " (meta)" : ""));
-      
-        Checkpoint ("Scheduling CPINITIAL");
-        CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction);
-      
-        Checkpoint ("Scheduling ANALYSIS");
-        CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
-      
-        Checkpoint ("OutputGH");
-        CCTK_OutputGH (cgh);
-      
-        PoisonCheck (cgh, alltimes);
-        CheckChecksums (cgh, allbutcurrenttime);
-      
+
+        const int do_every =
+          ipow(mgfact, mglevel) * (maxreflevelfact / ipow(reffact, rl));
+        if (cgh->cctk_iteration % do_every == 0)
+        {
+          // Checkpoint
+          Checkpoint ("Scheduling CPINITIAL");
+          CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction);
+
+          // Analysis
+          Checkpoint ("Scheduling ANALYSIS");
+          CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction);
+
+          // Output
+          Checkpoint ("OutputGH");
+          CCTK_OutputGH (cgh);
+
+          // Checking
+          PoisonCheck (cgh, alltimes);
+          CheckChecksums (cgh, allbutcurrenttime);
+        }
+
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
     }