Add Subcycling Feature

CarpetX/param.ccl

@@ -82,6 +82,13 @@ BOOLEAN periodic_z "Periodic"
 } no
 
 
+RESTRICTED:
+
+BOOLEAN use_subcycling_wip "Use subcycling in time"
+{
+} no
+
+
 
 PRIVATE:
 
@@ -444,12 +451,6 @@ CCTK_INT interpolation_order "Interpolation order" STEERABLE=always
 
 
 
-BOOLEAN use_subcycling_wip "Use subcycling in time"
-{
-} no  
-
-
-
 BOOLEAN do_reflux "Manage flux registers to ensure conservation"
 {
 } yes

CarpetX/src/driver.cxx

@@ -1777,20 +1777,29 @@ void CactusAmrCore::MakeNewLevelFromCoarse(
 #pragma omp critical
     CCTK_VINFO("MakeNewLevelFromCoarse patch %d level %d", patch, level);
 
-  assert(!use_subcycling_wip);
   assert(level > 0);
 
   SetupLevel(level, ba, dm, []() { return "MakeNewLevelFromCoarse"; });
 
   // Prolongate
-  assert(!use_subcycling_wip);
   auto &patchdata = ghext->patchdata.at(patch);
   auto &leveldata = patchdata.leveldata.at(level);
   auto &coarseleveldata = patchdata.leveldata.at(level - 1);
   const active_levels_t active_levels(level, level + 1, patch, patch + 1);
   const active_levels_t active_coarse_levels(level - 1, level, patch,
                                              patch + 1);
 
+  // only allow creation of new levels when the source and target are aligned
+  // in time
+  if (leveldata.iteration != coarseleveldata.iteration) {
+    ostringstream msg;
+    msg << "Coarse (rl=" << (level-1) <<", it=" << coarseleveldata.iteration <<
+           ") and fine (rl=" << level << ", it=" << leveldata.iteration <<
+           ") grid do not align in time when regridding";
+#pragma omp critical
+    CCTK_VERROR(msg.str().c_str());
+  }
+
   const int num_groups = CCTK_NumGroups();
   for (int gi = 0; gi < num_groups; ++gi) {
     cGroup group;
@@ -1886,8 +1895,16 @@ void CactusAmrCore::RemakeLevel(const int level, const amrex::Real time,
   const active_levels_t active_coarse_levels(level - 1, level, patch,
                                              patch + 1);
 
-  // Copy or prolongate
-  assert(!use_subcycling_wip);
+  // only allow modfication of levels when the source and target are aligned in
+  // time
+  if (leveldata.iteration != coarseleveldata.iteration) {
+    ostringstream msg;
+    msg << "Coarse (rl=" << (level-1) <<", it=" << coarseleveldata.iteration <<
+           ") and fine (rl=" << level << ", it=" << leveldata.iteration <<
+           ") grid do not align in time when regridding";
+#pragma omp critical
+    CCTK_VERROR(msg.str().c_str());
+  }
 
   // Check old level
   const int num_groups = CCTK_NumGroups();

CarpetX/src/schedule.cxx

@@ -83,6 +83,7 @@ std::optional<active_levels_t> active_levels;
 void Reflux(const cGH *cctkGH, int level);
 void Restrict(const cGH *cctkGH, int level, const std::vector<int> &groups);
 void Restrict(const cGH *cctkGH, int level);
+void SyncAfterRestrict(const cGH *cctkGH);
 
 namespace {
 // Convert a (direction, face) pair to an AMReX Orientation
@@ -492,6 +493,7 @@ void update_cctkGH(cGH *const cctkGH, const cGH *const sourceGH) {
   if (cctkGH == sourceGH)
     return;
   cctkGH->cctk_iteration = sourceGH->cctk_iteration;
+  cctkGH->cctk_timefac = sourceGH->cctk_timefac;
   cctkGH->cctk_time = sourceGH->cctk_time;
   cctkGH->cctk_delta_time = sourceGH->cctk_delta_time;
   // for (int d = 0; d < dim; ++d)
@@ -1352,6 +1354,8 @@ int Initialise(tFleshConfig *config) {
       if (leveldata.level != ghext->num_levels() - 1)
         Restrict(cctkGH, leveldata.level);
     });
+    // Prolongation
+    SyncAfterRestrict(cctkGH);
     CCTK_Traverse(cctkGH, "CCTK_POSTRESTRICT");
   }
 
@@ -1473,8 +1477,6 @@ void CycleTimelevels(cGH *restrict const cctkGH) {
   static Timer timer("CycleTimelevels");
   Interval interval(timer);
 
-  cctkGH->cctk_iteration += 1;
-  cctkGH->cctk_time += cctkGH->cctk_delta_time;
   update_cctkGHs(cctkGH);
 
   // TODO: Parallelize over groups
@@ -1712,13 +1714,14 @@ int Evolve(tFleshConfig *config) {
       for (const auto &leveldata : patchdata.leveldata)
         iteration = min(iteration, leveldata.iteration);
 
+    cctkGH->cctk_iteration += 1;
+
     // Loop over all levels, in batches that combine levels that don't
     // subcycle. The level range is [min_level, max_level).
-    int min_level = 0;
-    while (min_level < ghext->num_levels()) {
+    for (int min_level = 0, max_level = min_level + 1;
+         min_level < ghext->num_levels();
+         min_level = max_level, max_level = min_level + 1) {
       // Find end of batch
-      int max_level = min_level + 1;
-
       while (max_level < ghext->num_levels()) {
         bool level_is_subcycling_level = false;
         for (const auto &patchdata : ghext->patchdata)
@@ -1733,18 +1736,31 @@ int Evolve(tFleshConfig *config) {
       // Skip this batch of levels if it is not active at the current
       // iteration
       rat64 level_iteration = -1;
-      for (const auto &patchdata : ghext->patchdata)
-        if (min_level < int(patchdata.leveldata.size()))
+      rat64 level_delta_iteration = -1;
+      for (const auto &patchdata : ghext->patchdata) {
+        if (min_level < int(patchdata.leveldata.size())) {
           level_iteration = patchdata.leveldata.at(min_level).iteration;
+          level_delta_iteration =
+              patchdata.leveldata.at(min_level).delta_iteration;
+          break;
+        }
+      }
       assert(level_iteration != -1);
+      assert(level_delta_iteration != -1);
+      // Skip those evolved coarse levels when evolving fine levels to catch up
       if (level_iteration > iteration)
-        break;
+        continue;
 
+      // must not terminate loop iteration due to active_levels being reset at
+      // bootom of loop body
       active_levels = make_optional<active_levels_t>(min_level, max_level);
 
       // Advance iteration number on this batch of levels
+      level_iteration += level_delta_iteration;
       active_levels->loop_serially([&](auto &restrict leveldata) {
         leveldata.iteration += leveldata.delta_iteration;
+        assert(level_iteration == leveldata.iteration);
+        assert(level_delta_iteration == leveldata.delta_iteration);
       });
 
       // We cannot invalidate all non-evolved variables. ODESolvers
@@ -1754,6 +1770,11 @@ int Evolve(tFleshConfig *config) {
 
       CycleTimelevels(cctkGH);
 
+      cctkGH->cctk_timefac = use_subcycling_wip ? std::pow(2, min_level) : 1;
+      cctkGH->cctk_time =
+          use_subcycling_wip ? cctkGH->cctk_delta_time * double(level_iteration)
+                             : cctkGH->cctk_time + cctkGH->cctk_delta_time;
+
       CCTK_Traverse(cctkGH, "CCTK_PRESTEP");
       CCTK_Traverse(cctkGH, "CCTK_EVOL");
 
@@ -1763,11 +1784,25 @@ int Evolve(tFleshConfig *config) {
       for (int level = ghext->num_levels() - 2; level >= 0; --level)
         Reflux(cctkGH, level);
 
+      // reset active_levels to all that have caught to the same time
+      for (int level = min_level - 1; level >= 0; --level) {
+        rat64 coarse_iteration =
+            ghext->patchdata.at(0).leveldata.at(level).iteration;
+        if (coarse_iteration == level_iteration)
+          min_level = level;
+        else
+          break;
+      }
+      active_levels = make_optional<active_levels_t>(min_level, max_level);
+
       if (!restrict_during_sync) {
         // Restrict
-        // TODO: These loop bounds are wrong for subcycling
-        for (int level = ghext->num_levels() - 2; level >= 0; --level)
-          Restrict(cctkGH, level);
+        active_levels->loop_fine_to_coarse([&](const auto &leveldata) {
+          if (leveldata.level + 1 < active_levels->max_level)
+            Restrict(cctkGH, leveldata.level);
+        });
+        // Prolongation
+        SyncAfterRestrict(cctkGH);
         CCTK_Traverse(cctkGH, "CCTK_POSTRESTRICT");
       }
 
@@ -1780,7 +1815,7 @@ int Evolve(tFleshConfig *config) {
       total_evolution_output_time += output_finish_time - output_start_time;
 
       active_levels = optional<active_levels_t>();
-    } // for min_level
+    } // for min_level, max_level
 
     const double finish_time = gettime();
     double num_cells = 0;
@@ -2784,6 +2819,24 @@ void Restrict(const cGH *cctkGH, int level) {
   Restrict(cctkGH, level, groups);
 }
 
+void SyncAfterRestrict(const cGH *cctkGH) {
+  const int numgroups = CCTK_NumGroups();
+  vector<int> groups;
+  groups.reserve(numgroups);
+  const auto &patchdata0 = ghext->patchdata.at(0);
+  const auto &leveldata0 = patchdata0.leveldata.at(0);
+  for (const auto &groupdataptr : leveldata0.groupdata) {
+    // Sync only grid functions
+    if (groupdataptr) {
+      auto &restrict groupdata = *groupdataptr;
+      // Sync only evolved grid functions
+      if (groupdata.do_checkpoint)
+        groups.push_back(groupdata.groupindex);
+    }
+  }
+  SyncGroupsByDirI(cctkGH, groups.size(), groups.data(), nullptr);
+}
+
 // storage handling
 namespace {
 int GroupStorageCrease(const cGH *cctkGH, int n_groups, const int *groups,

TestSubcycling/README

@@ -0,0 +1,9 @@
+Cactus Code Thorn TestSubcycling
+Author(s)    : Roland Haas <[email protected]>
+Maintainer(s): Roland Haas <[email protected]>
+Licence      : NCSA
+--------------------------------------------------------------------------
+
+1. Purpose
+
+not documented

TestSubcycling/configuration.ccl

@@ -0,0 +1,3 @@
+# Configuration definitions for thorn TestSubcycling
+
+REQUIRES Loop