Carpet: Store number of active time levels per level instead of globally

Store the number of active time levels explicitly for each refinement
level and convergence level.

darcs-hash:20060413192937-dae7b-579ddb7ae741d73bb8bf6e0c87a4bff1c18f7c46.gz

6 files changed, 279 insertions, 80 deletions

diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index c1e01cf9c..063dd20f1 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -165,6 +165,18 @@ namespace Carpet {
 
     } END_MAP_LOOP;
     
+    // Set new number of active time levels
+    for (int n=0; n<CCTK_NumGroups(); ++n) {
+      int const grouptype = CCTK_GroupTypeI (n);
+      if (grouptype == CCTK_GF) {
+        for (int ml=0; ml<mglevels; ++ml) {
+          groupdata.at(n).activetimelevels.at(ml).resize
+            (vhh.at(0)->reflevels(),
+             groupdata.at(n).activetimelevels.at(ml).at(0));
+        }
+      }
+    }
+    
     // Calculate new number of levels
     reflevels = vhh.at(0)->reflevels();
     for (int m=0; m<maps; ++m) {
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 0ef292fc3..a70b89ebb 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -240,14 +240,24 @@ namespace Carpet {
       if (can_transfer) {
         // Use the default
         if (gp.numtimelevels == 1) {
+#if 0
           // Only one time level: copy instead of prolongating
           char * const groupname = CCTK_GroupName (group);
           CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                       "Group \"%s\" has only one time level; therefore it "
-                      "will be copied insted of prolongated.",
+                      "will be copied instead of prolongated.",
                       groupname);
           free (groupname);
           return op_copy;
+#endif
+          // Only one time level:
+          char * const groupname = CCTK_GroupName (group);
+          CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "Group \"%s\" has only one time level; therefore it "
+                      "will not be prolongated or restricted.",
+                      groupname);
+          free (groupname);
+          return op_none;
         } else {
           // Several time levels: use the default
           return op_Lagrange;
@@ -1133,6 +1143,10 @@ namespace Carpet {
 
     case CCTK_GF: {
       assert (gp.dim == dim);
+      groupdata.at(group).activetimelevels.resize(mglevels);
+      for (int ml=0; ml<mglevels; ++ml) {
+        groupdata.at(group).activetimelevels.at(ml).resize(1);
+      }
       arrdata.at(group).resize(maps);
       for (int m=0; m<maps; ++m) {
         arrdata.at(group).at(m).hh = vhh.at(m);
@@ -1145,6 +1159,10 @@ namespace Carpet {
     case CCTK_SCALAR:
     case CCTK_ARRAY: {
         
+      groupdata.at(group).activetimelevels.resize(mglevels);
+      for (int ml=0; ml<mglevels; ++ml) {
+        groupdata.at(group).activetimelevels.at(ml).resize(1);
+      }
       arrdata.at(group).resize(1);
       
       ivect sizes(1), ghostsizes(0);
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index a48db2a7a..257154704 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -56,6 +56,14 @@ namespace Carpet {
     for (int n=0; n<n_groups; ++n) {
       int const group = groups[n];
       
+      {
+        char * const groupname = CCTK_GroupName (group);
+        assert (groupname);
+        Checkpoint ("  %s: %screase to %d",
+                    groupname, inc ? "in" : "de", timelevels[n]);
+        free (groupname);
+      }
+      
       cGroup gp;
       const int ierr = CCTK_GroupData (group, &gp);
       assert (not ierr);
@@ -72,49 +80,51 @@ namespace Carpet {
       // Check an assumption
       if (not gp.vectorgroup) assert (gp.vectorlength == 1);
       
-      // Record previous number of allocated time levels
-      if (status) {
-        status[n] = groupdata.at(group).info.activetimelevels;
+      // No storage change in local mode
+      if (gp.grouptype == CCTK_GF) {
+        assert ((map == -1 or maps == 1)
+                and (component == -1
+                     or vhh.at(0)->local_components(reflevel) == 1));
       }
       
-      // Only do something if the number of time levels actually needs
-      // to be changed -- do nothing otherwise
-
-      const bool do_increase
-        = inc and timelevels[n] > groupdata.at(group).info.activetimelevels;
-      const bool do_decrease
-        = ! inc and timelevels[n] < groupdata.at(group).info.activetimelevels;
-      if (do_increase or do_decrease) {
-        
-        // No storage change in local mode
-        if (gp.grouptype == CCTK_GF) {
-          assert ((map == -1 or maps == 1)
-                  and (component == -1
-                       or vhh.at(0)->local_components(reflevel) == 1));
-        }
-        
-        if (! can_do) {
-          char * const groupname = CCTK_GroupName (group);
-          char const * const modestring
-            = (is_meta_mode() ? "meta"
-               : is_global_mode() ? "global"
-               : is_level_mode() ? "level"
-               : is_singlemap_mode() ? "singlemap"
-               : is_local_mode() ? "local"
-               : NULL);
-          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                      "Cannot change storage for group \"%s\" in %s mode",
-                      groupname, modestring);
-          free (groupname);
-        }
-        assert (can_do);
-        
-        // Set the new number of active time levels
-        groupdata.at(group).info.activetimelevels = timelevels[n];
-        
-        // Allocate the time levels as well
-        for (int ml=min_ml; ml<max_ml; ++ml) {
-          for (int rl=min_rl; rl<max_rl; ++rl) {
+      // Allocate the time levels
+      for (int ml=min_ml; ml<max_ml; ++ml) {
+        for (int rl=min_rl; rl<max_rl; ++rl) {
+          
+          // Record previous number of allocated time levels
+          if (status) {
+            // Note: This remembers only the last level
+            status[n] = groupdata.at(group).activetimelevels.at(ml).at(rl);
+          }
+          
+          // Only do something if the number of time levels actually
+          // needs to be changed -- do nothing otherwise
+          
+          const bool do_increase
+            =     inc and timelevels[n] > groupdata.at(group).activetimelevels.at(ml).at(rl);
+          const bool do_decrease
+            = not inc and timelevels[n] < groupdata.at(group).activetimelevels.at(ml).at(rl);
+          if (do_increase or do_decrease) {
+            
+            if (not can_do) {
+              char * const groupname = CCTK_GroupName (group);
+              char const * const modestring
+                = (is_meta_mode() ? "meta"
+                   : is_global_mode() ? "global"
+                   : is_level_mode() ? "level"
+                   : is_singlemap_mode() ? "singlemap"
+                   : is_local_mode() ? "local"
+                   : NULL);
+              CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                          "Cannot change storage for group \"%s\" in %s mode",
+                          groupname, modestring);
+              free (groupname);
+            }
+            assert (can_do);
+            
+            // Set the new number of active time levels
+            groupdata.at(group).activetimelevels.at(ml).at(rl) = timelevels[n];
+            
             for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
               for (int var=0; var<gp.numvars; ++var) {
                 const int vectorindex
@@ -157,7 +167,7 @@ namespace Carpet {
                   int const c = CCTK_MyProc(cgh);
                   for (int tl=0; tl<gp.numtimelevels; ++tl) {
                     cgh->data[varindex][tl]
-                      = (tl < groupdata.at(group).info.activetimelevels
+                      = (tl < groupdata.at(group).activetimelevels.at(ml).at(rl)
                          ? ((*arrdata.at(group).at(m).data.at(var))
                             (tl, 0, c, 0)->storage())
                          : NULL);
@@ -166,42 +176,45 @@ namespace Carpet {
               
               } // for var
             } // for m
-          } // for rl
-        } // for ml
-        
-      } // if really change the number of active time levels
-      
-      // Complain if there are not enough active time levels
-      if (gp.grouptype == CCTK_GF) {
-        if (max_refinement_levels > 1) {
-          if (groupdata.at(group).transport_operator != op_none
-              and groupdata.at(group).transport_operator != op_copy) {
-            if (groupdata.at(group).info.activetimelevels != 0
-                and (groupdata.at(group).info.activetimelevels
-                     < prolongation_order_time+1))
-            {
-              static vector<bool> didwarn;
-              int const numgroups = CCTK_NumGroups();
-              if (didwarn.size() < numgroups) {
-                didwarn.resize (numgroups, false);
-              }
-              if (! didwarn.at(group)) {
-                // Warn only once per group
-                didwarn.at(group) = true;
-                char * const groupname = CCTK_GroupName (group);
-                CCTK_VWarn (1, __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.",
-                            groupname,
-                            prolongation_order_time, prolongation_order_time+1);
-                free (groupname);
+            
+          } // if really change the number of active time levels
+          
+          // Complain if there are not enough active time levels
+          if (gp.grouptype == CCTK_GF) {
+            if (max_refinement_levels > 1) {
+              if (groupdata.at(group).transport_operator != op_none
+                  and groupdata.at(group).transport_operator != op_copy) {
+                if (groupdata.at(group).activetimelevels.at(ml).at(rl) != 0
+                    and (groupdata.at(group).activetimelevels.at(ml).at(rl)
+                         < prolongation_order_time+1))
+                {
+                  static vector<bool> didwarn;
+                  int const numgroups = CCTK_NumGroups();
+                  if (didwarn.size() < numgroups) {
+                    didwarn.resize (numgroups, false);
+                  }
+                  if (not didwarn.at(group)) {
+                    // Warn only once per group
+                    didwarn.at(group) = true;
+                    char * const groupname = CCTK_GroupName (group);
+                    CCTK_VWarn (1, __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.",
+                                groupname,
+                                prolongation_order_time, prolongation_order_time+1);
+                    free (groupname);
+                  }
+                }
               }
             }
           }
-        }
-      }
-      
-      // Record current number of time levels
-      total_num_timelevels += groupdata.at(group).info.activetimelevels;
+          
+          // Record current number of time levels
+          // Note: This adds the time levels of all refinement levels
+          total_num_timelevels
+            += groupdata.at(group).activetimelevels.at(ml).at(rl);
+          
+        } // for rl
+      } // for ml
       
     } // for n
     
@@ -263,12 +276,17 @@ namespace Carpet {
   int
   QueryGroupStorageB (const cGH* cgh, int group, const char* groupname)
   {
+    DECLARE_CCTK_PARAMETERS;
     if (groupname) {
       group = CCTK_GroupIndex(groupname);
     }
     assert (group>=0 and group<CCTK_NumGroups());
+    int const grouptype = CCTK_GroupTypeI (group);
+    if (is_meta_mode() or is_global_mode()) assert (grouptype != CCTK_GF);
+    int const rl = grouptype == CCTK_GF ? reflevel : 0;
     // Return whether storage is allocated
-    return groupdata.at(group).info.activetimelevels > 0;
+    assert (groupdata.at(group).activetimelevels.at(mglevel).at(rl) != deadbeef);
+    return groupdata.at(group).activetimelevels.at(mglevel).at(rl) > 0;
   }
   
   
diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc
index bd674ab83..bb0df526e 100644
--- a/Carpet/Carpet/src/modes.cc
+++ b/Carpet/Carpet/src/modes.cc
@@ -111,6 +111,8 @@ namespace Carpet {
           const_cast<int*>(groupdata.at(group).info.bbox)[2*d  ] = obnds[d][0];
           const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1] = obnds[d][1];
         }
+        groupdata.at(group).info.activetimelevels
+          = groupdata.at(group).activetimelevels.at(mglevel).at(0);
         
         for (int d=0; d<dim; ++d) {
           assert (groupdata.at(group).info.lsh[d]>=0);
@@ -129,7 +131,7 @@ namespace Carpet {
           assert (firstvar>=0);
           const int max_tl = CCTK_MaxTimeLevelsGI (group);
           assert (max_tl>=0);
-          const int active_tl = CCTK_ActiveTimeLevelsGI (cgh, group);
+          const int active_tl = groupdata.at(group).info.activetimelevels;
           assert (active_tl>=0 and active_tl<=max_tl);
           
           assert (arrdata.at(group).at(m).hh->is_local(rl,c));
@@ -199,6 +201,7 @@ namespace Carpet {
           const_cast<int*>(groupdata.at(group).info.bbox)[2*d  ] = deadbeef;
           const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1] = deadbeef;
         }
+        groupdata.at(group).info.activetimelevels = deadbeef;
         
         const int numvars = CCTK_NumVarsInGroupI (group);
         if (numvars>0) {
@@ -244,6 +247,14 @@ namespace Carpet {
     ivect::ref(cgh->cctk_levfac) = spacereflevelfact;
     cgh->cctk_timefac = timereflevelfact;
     
+    // Set number of time levels
+    for (int group=0; group<CCTK_NumGroups(); ++group) {
+      if (CCTK_GroupTypeI(group) == CCTK_GF) {
+        groupdata.at(group).info.activetimelevels
+          = groupdata.at(group).activetimelevels.at(mglevel).at(reflevel);
+      }
+    }
+    
     // Set current time
     assert (mglevel>=0 and mglevel<(int)leveltimes.size());
     assert (reflevel>=0 and reflevel<(int)leveltimes.at(mglevel).size());
@@ -275,6 +286,13 @@ namespace Carpet {
       global_time = cgh->cctk_time;
     }
     
+    // Unset number of time levels
+    for (int group=0; group<CCTK_NumGroups(); ++group) {
+      if (CCTK_GroupTypeI(group) == CCTK_GF) {
+        groupdata.at(group).info.activetimelevels = deadbeef;
+      }
+    }
+    
     reflevel = -1;
     timereflevelfact = timereffacts.at (reflevels - 1);
     // TODO: use spacereffacts.at (reflevel - 1) instead?
@@ -1008,6 +1026,68 @@ namespace Carpet {
   
   
   //
+  // Mode setters
+  //
+  
+  // mglevel setter
+  
+  mglevel_setter::mglevel_setter (cGH const * const cgh_, int const ml)
+    : cgh(const_cast<cGH*>(cgh_))
+  {
+    assert (is_meta_mode());
+    enter_global_mode (cgh, ml);
+  }
+  
+  mglevel_setter::~mglevel_setter ()
+  {
+    leave_global_mode (cgh);
+  }
+  
+  // reflevel setter
+  
+  reflevel_setter::reflevel_setter (cGH const * const cgh_, int const rl)
+    : cgh(const_cast<cGH*>(cgh_))
+  {
+    assert (is_global_mode());
+    enter_level_mode (cgh, rl);
+  }
+  
+  reflevel_setter::~reflevel_setter ()
+  {
+    leave_level_mode (cgh);
+  }
+  
+  // map setter
+  
+  map_setter::map_setter (cGH const * const cgh_, int const m)
+    : cgh(const_cast<cGH*>(cgh_))
+  {
+    assert (is_level_mode());
+    enter_singlemap_mode (cgh, m);
+  }
+  
+  map_setter::~map_setter ()
+  {
+    leave_singlemap_mode (cgh);
+  }
+  
+  // component setter
+  
+  component_setter::component_setter (cGH const * const cgh_, int const c)
+    : cgh(const_cast<cGH*>(cgh_))
+  {
+    assert (is_singlemap_mode());
+    enter_local_mode (cgh, c);
+  }
+  
+  component_setter::~component_setter ()
+  {
+    leave_local_mode (cgh);
+  }
+  
+  
+  
+  //
   // Call a scheduling group
   //
   
diff --git a/Carpet/Carpet/src/modes.hh b/Carpet/Carpet/src/modes.hh
index 001e1cc7e..3b1cd873e 100644
--- a/Carpet/Carpet/src/modes.hh
+++ b/Carpet/Carpet/src/modes.hh
@@ -136,8 +136,8 @@ namespace Carpet {
     bool done () const;
     void step ();
   };
-    
-    
+  
+  
   
   // Compatibility defines for the mode iterators
   
@@ -315,6 +315,76 @@ namespace Carpet {
     meta_mode_ = false;                         \
   } while (false)
   
+  
+  
+  // Mode setters
+  
+  class mglevel_setter {
+    cGH * cgh;
+  public:
+    mglevel_setter (cGH const * const cgh, int const ml);
+    ~mglevel_setter ();
+  };
+  
+  class reflevel_setter {
+    cGH * cgh;
+  public:
+    reflevel_setter (cGH const * const cgh, int const rl);
+    ~reflevel_setter ();
+  };
+  
+  class map_setter {
+    cGH * cgh;
+  public:
+    map_setter (cGH const * const cgh, int const m);
+    ~map_setter ();
+  };
+  
+  class component_setter {
+    cGH * cgh;
+  public:
+    component_setter (cGH const * const cgh, int const c);
+    ~component_setter ();
+  };
+  
+  
+  
+  // Compatibility defines for the mode setters
+  
+#define ENTER_GLOBAL_MODE(cgh, ml)              \
+  do {                                          \
+    mglevel_setter mg_setter_(cgh, ml);         \
+    {
+#define LEAVE_GLOBAL_MODE                       \
+    }                                           \
+  } while (false)
+  
+#define ENTER_LEVEL_MODE(cgh, rl)               \
+  do {                                          \
+    reflevel_setter ref_setter_(cgh, rl);       \
+    {
+#define LEAVE_LEVEL_MODE                        \
+    }                                           \
+  } while (false)
+  
+#define ENTER_SINGLEMAP_MODE(cgh, m)            \
+  do {                                          \
+    map_setter m_setter_(cgh, m);               \
+    {
+#define LEAVE_SINGLEMAP_MODE                    \
+    }                                           \
+  } while (false)
+  
+#define ENTER_LOCAL_MODE(cgh, c)                \
+  do {                                          \
+    component_setter c_setter_(cgh, c);         \
+    {
+#define LEAVE_LOCAL_MODE                        \
+    }                                           \
+  } while (false)
+  
+  
+  
 } // namespace Carpet
 
 #endif // !defined(MODES_HH)
diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh
index b2981ee24..94bb57843 100644
--- a/Carpet/Carpet/src/variables.hh
+++ b/Carpet/Carpet/src/variables.hh
@@ -129,6 +129,7 @@ namespace Carpet {
   struct groupdesc {
     cGroupDynamicData info;
     operator_type transport_operator; // prolongation and restriction
+    vector<vector<int> > activetimelevels; // [mglevel][reflevel]
   };
   extern vector<groupdesc> groupdata; // [group]