Add more global time step control

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1159]](https://github.com/parthenon-hpc-lab/parthenon/pull/1159) Add additional timestep controllers in parthenon/time.
 - [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
 - [[PR 1143]](https://github.com/parthenon-hpc-lab/parthenon/pull/1143) Add tensor indices to VariableState, add radiation constant to constants, add TypeLists, allow for arbitrary containers for solvers
 - [[PR 1140]](https://github.com/parthenon-hpc-lab/parthenon/pull/1140) Allow for relative convergence tolerance in BiCGSTAB solver.

doc/sphinx/src/inputs.rst

@@ -21,8 +21,8 @@ General parthenon options such as problem name and parameter handling.
 +---------------------+---------+---------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Option              | Default | Type    | Description                                                                                                                                                                                            |
 +=====================+=========+=========+========================================================================================================================================================================================================+
-|| name               || none   || string || Name of this problem or initialization, prefixed to output files.                                                                                                                                     |
 || archive_parameters || false  || string || Produce a parameter file containing all parameters known to Parthenon. Set to `true` for an output file named `parthinput.archive`. Set to `timestamp` for a file with a name containing a timestamp. |
+|| name               || none   || string || Name of this problem or initialization, prefixed to output files.                                                                                                                                     |
 +---------------------+---------+---------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 
 
@@ -34,12 +34,23 @@ Options related to time-stepping and printing of diagnostic data.
 +------------------------------+---------+--------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Option                       | Default | Type   | Description                                                                                                                                                           |
 +==============================+=========+========+=======================================================================================================================================================================+
-|| tlim                        || none   || float || Stop criterion on simulation time.                                                                                                                                   |
-|| nlim                        || -1     || int   || Stop criterion on total number of steps taken. Ignored if < 0.                                                                                                       |
-|| perf_cycle_offset           || 0      || int   || Skip the first N cycles when calculating the final performance (e.g., zone-cycles/wall_second). Allows to hide the initialization overhead in Parthenon.             |
+|| dt_ceil                     || none   || Real  || The maximum allowed timestep.                                                                                                                                        |
+|| dt_factor                   || 2.0    || Real  || The maximum allowed relative increase of the timestep over the previous value.                                                                                       |
+|| dt_floor                    || none   || Real  || The minimum allowed timestep.                                                                                                                                        |
+|| dt_force                    || none   || Real  || Force the timestep to this value, ignoring all other conditions.                                                                                                     |
+|| dt_init                     || none   || Real  || The maximum allowed timestep during the first cycle.                                                                                                                 |
+|| dt_init_force               || none   || bool  || If set to true, force the first cycle's timestep to the value given by dt_init.                                                                                      |
+|| dt_min                      || none   || Real  || If the timestep falls below dt_min for dt_min_cycle_limit cycles, Parthenon fatals.                                                                                  |
+|| dt_min_cycle_limit          || 10     || int   || The maximum number of cycles the timestep can be below dt_min.                                                                                                       |
+|| dt_max                      || none   || Real  || If the timestep falls below dt_max for dt_max_cycle_limit cycles, Parthenon fatals.                                                                                  |
+|| dt_max_cycle_limit          || 1      || int   || The maximum number of cycles the timestep an be above dt_max.                                                                                                        |
+|| dt_user                     || none   || Real  || Set a global timestep limit.                                                                                                                                         |
+|| ncrecv_bdry_buf_timeout_sec || -1.0   || Real  || Timeout in seconds for the `ReceiveBoundaryBuffers` tasks. Disabed (negative) by default. Typically no need in production runs. Useful for debugging MPI calls.      |
 || ncycle_out                  || 1      || int   || Number of cycles between short diagnostic output to standard out containing, e.g., current time, dt, zone-update/wsec. Default: 1 (i.e, every cycle).                |
 || ncycle_out_mesh             || 0      || int   || Number of cycles between printing the mesh structure to standard out. Use a negative number to also print every time the mesh was modified. Default: 0 (i.e, off).   |
-|| ncrecv_bdry_buf_timeout_sec || -1.0   || Real  || Timeout in seconds for the `ReceiveBoundaryBuffers` tasks. Disabed (negative) by default. Typically no need in production runs. Useful for debugging MPI calls.      |
+|| nlim                        || -1     || int   || Stop criterion on total number of steps taken. Ignored if < 0.                                                                                                       |
+|| perf_cycle_offset           || 0      || int   || Skip the first N cycles when calculating the final performance (e.g., zone-cycles/wall_second). Allows to hide the initialization overhead in Parthenon.             |
+|| tlim                        || none   || Real  || Stop criterion on simulation time.                                                                                                                                   |
 +------------------------------+---------+--------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 
 
@@ -64,9 +75,9 @@ See the :ref:`sparse impl` documentation for details.
 +--------------------+---------+--------+----------------------------------------------------------------------------------------------------------------------------------------------+
 | Option             | Default | Type   | Description                                                                                                                                  |
 +====================+=========+========+==============================================================================================================================================+
-|| enable_sparse     || `true` || bool  || If set to false, sparse variables will always be allocated, see also :ref:`sparse run-time`                                                 |
 || alloc_threshold   || 1e-12  || float || Global (for all sparse variables) threshold to trigger allocation of a variable if cells in the receiving ghost cells are above this value. |
-|| dealloc_threshold || 1e-14  || float || Global (for all sparse variables) threshold to trigger deallocation if all active cells of a variable in a block are below this value.      |
 || dealloc_count     || 5      || int   || First deallocate a sparse variable if the `dealloc_threshold` has been met in this number of consecutive cycles.                            |
+|| dealloc_threshold || 1e-14  || float || Global (for all sparse variables) threshold to trigger deallocation if all active cells of a variable in a block are below this value.      |
+|| enable_sparse     || `true` || bool  || If set to false, sparse variables will always be allocated, see also :ref:`sparse run-time`                                                 |
 +--------------------+---------+--------+----------------------------------------------------------------------------------------------------------------------------------------------+
 

src/driver/driver.cpp

@@ -205,25 +205,65 @@ void EvolutionDriver::InitializeBlockTimeSteps() {
 // \brief function that loops over all MeshBlocks and find new timestep
 
 void EvolutionDriver::SetGlobalTimeStep() {
-  // don't allow dt to grow by more than 2x
-  // consider making this configurable in the input
-  if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
-    tm.dt *= 2.0;
-  }
-  Real big = std::numeric_limits<Real>::max();
-  for (auto const &pmb : pmesh->block_list) {
-    tm.dt = std::min(tm.dt, pmb->NewDt());
-    pmb->SetAllowedDt(big);
-  }
-
+  if (dt_force > 0.0) {
+    // Check if user wants to force the value
+    tm.dt = dt_force;
+  } else if ((tm.ncycle == 0) && (dt_init_force) && (dt_init > 0.0)) {
+    // Check if user wants to force the first cycle value
+    tm.dt = dt_init;
+  } else {
+    if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
+      // don't allow dt to grow by more than 2x
+      // consider making this configurable in the input
+      tm.dt *= dt_factor;
+    }
+    // Check for special first cycle value
+    if (tm.ncycle == 0) {
+      tm.dt = std::min(tm.dt, dt_init);
+    }
+    // Allow the meshblocks to vote
+    for (auto const &pmb : pmesh->block_list) {
+      tm.dt = std::min(tm.dt, pmb->NewDt());
+      pmb->SetAllowedDt(std::numeric_limits<Real>::max());
+    }
+    // Allow the user to enforce maximum timestep
+    tm.dt = std::min(tm.dt, dt_user);
+    // Force timestep to be in the allowable range
+    tm.dt = std::max(dt_floor, std::min(tm.dt, dt_ceil));
 #ifdef MPI_PARALLEL
-  PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &tm.dt, 1, MPI_PARTHENON_REAL, MPI_MIN,
-                                    MPI_COMM_WORLD));
+    PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &tm.dt, 1, MPI_PARTHENON_REAL,
+                                      MPI_MIN, MPI_COMM_WORLD));
 #endif
+  }
 
-  if (tm.time < tm.tlim &&
-      (tm.tlim - tm.time) < tm.dt) // timestep would take us past desired endpoint
+  // Check that we have not gone off the rails
+  if (tm.dt <= dt_min) {
+    if (++dt_min_count >= dt_min_count_max) {
+      std::stringstream msg;
+      msg << "Timesetep has fallen bellow minimum (parthenon/time/dt_min=" << dt_min
+          << ") for more than " << dt_min_count_max << " steps";
+      PARTHENON_FAIL(msg);
+    }
+  } else {
+    dt_min_count = 0;
+  }
+  if (tm.dt >= dt_max) {
+    if (++dt_max_count >= dt_max_count_max) {
+      std::stringstream msg;
+      msg << "Timesetep has risen above maximum (parthenon/time/dt_max=" << dt_max
+          << ") for more than " << dt_max_count_max << " steps";
+      PARTHENON_FAIL(msg);
+    }
+  } else {
+    dt_max_count = 0;
+  }
+
+  // Limit timestep if it would take us past desired endpoint
+  // This comes after the bounds checking so that we don't fail when epsilon
+  // away from tlim
+  if (tm.time < tm.tlim && (tm.tlim - tm.time) < tm.dt) {
     tm.dt = tm.tlim - tm.time;
+  }
 }
 
 void EvolutionDriver::DumpInputParameters() {

src/driver/driver.hpp

@@ -69,6 +69,30 @@ class EvolutionDriver : public Driver {
                                       std::numeric_limits<Real>::infinity());
     Real dt =
         pinput->GetOrAddReal("parthenon/time", "dt", std::numeric_limits<Real>::max());
+    dt_min = pinput->GetOrAddReal("parthenon/time", "dt_min",
+                                  std::numeric_limits<Real>::min());
+    dt_max = pinput->GetOrAddReal("parthenon/time", "dt_max",
+                                  std::numeric_limits<Real>::max());
+    dt_init = pinput->GetOrAddReal("parthenon/time", "dt_init",
+                                   std::numeric_limits<Real>::max());
+    dt_init_force = pinput->GetOrAddBoolean("parthenon/time", "dt_init_force", false);
+
+    dt_user = pinput->GetOrAddReal("parthenon/time", "dt_user",
+                                   std::numeric_limits<Real>::max());
+    dt_force = pinput->GetOrAddReal("parthenon/time", "dt_force",
+                                    std::numeric_limits<Real>::lowest());
+    dt_floor = pinput->GetOrAddReal("parthenon/time", "dt_floor",
+                                    std::numeric_limits<Real>::min());
+    dt_ceil = pinput->GetOrAddReal("parthenon/time", "dt_ceil",
+                                   std::numeric_limits<Real>::max());
+    dt_min_count_max =
+        pinput->GetOrAddInteger("parthenon/time", "dt_min_cycle_limit", 10);
+    dt_max_count_max = pinput->GetOrAddInteger("parthenon/time", "dt_max_cycle_limit", 1);
+    dt_min_count = 0;
+    dt_max_count = 0;
+
+    dt_factor = pinput->GetOrAddReal("parthenon/time", "dt_factor", 2.0);
+
     const auto ncycle = pinput->GetOrAddInteger("parthenon/time", "ncycle", 0);
     const auto nmax = pinput->GetOrAddInteger("parthenon/time", "nlim", -1);
     const auto nout = pinput->GetOrAddInteger("parthenon/time", "ncycle_out", 1);
@@ -88,6 +112,11 @@ class EvolutionDriver : public Driver {
 
  protected:
   void PostExecute(DriverStatus status) override;
+  Real dt_user, dt_force, dt_init, dt_min, dt_max, dt_floor, dt_ceil;
+  Real dt_factor;
+  bool dt_init_force;
+  int dt_min_count, dt_max_count;
+  int dt_min_count_max, dt_max_count_max;
 
  private:
   void InitializeBlockTimeSteps();