Generalize Integrator Infrastructure

.gitignore

@@ -41,3 +41,7 @@ compile_commands.json
 # Files created by bad CMake hygiene
 a.out
 cmake_hdf5_test.o
+
+# Python artifacts
+*.pyc
+*.egg-info

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 840]](https://github.com/parthenon-hpc-lab/parthenon/pull/840) Generalized integrators infrastructure in a backwards compatible way
 - [[PR 810]](https://github.com/parthenon-hpc-lab/parthenon/pull/810) Add suport for Ascent in-situ visualization
 - [[PR 831]](https://github.com/parthenon-hpc-lab/parthenon/pull/831) Add set-based MetadataFlag logic
 - [[PR 803]](https://github.com/parthenon-hpc-lab/parthenon/pull/803) Add skeleton for sphinx docs

doc/sphinx/src/driver.rst

@@ -16,7 +16,7 @@ It has a single pure virtual member function called ``Execute`` that
 must be defined by a derived class and is intended to be called from
 ``main()``. A simple example of defining an application based on
 inheriting from this class can be found
-`here <../example/calculate_pi/pi_driver.hpp>`__.
+`here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop/example/calculate_pi/pi_driver.hpp>`__.
 
 EvolutionDriver
 ---------------
@@ -38,9 +38,10 @@ MultiStageDriver
 ----------------
 
 The ``MultiStageDriver`` derives from the ``EvolutionDriver``, extending
-it in two important ways. First, it holds a
-``std::unique_ptr<StagedIntegrator>`` object which includes members for
+it in two important ways. First, it is templated on an ``Integrator`` type and
+holds a ``std::unique_ptr<Integrator>`` object which includes members for
 the number of stages, the stage weights, and the names of the stages.
+More details on integrators can be found on the :ref:`integrators` page.
 Second, it defines a ``Step()`` function, which is reponsible for taking
 a timestep by looping over stages and calling the
 ``ConstructAndExecuteTaskLists`` function which builds and executes the

doc/sphinx/src/integrators.rst

@@ -0,0 +1,96 @@
+.. _integrators:
+
+Integrators
+============
+
+Integrators contain useful properties for writing a time-integration
+scheme. They are defined in a class hierarchy. ``MultiStageDriver``
+owns a pointer to an integrator. The type of integrator depends on the
+template parameter. By default it is the base class, and you can
+assign child integrator you like, but you must cast the pointer at
+runtime to access useful features of the child classes. Alternatively
+if you inherit from ``MultiStageDriver`` and specialize to a specific
+child class, you can avoid doing the cast.
+
+StagedIntegrator
+------------------
+
+The base class is the ``StagedIntegrator``. It contains a timestep
+``dt``, a total number of integration steps, ``nstages``, a total
+number of scratch buffers required, ``nbuffers``, and lists of strings
+containing suggested names for the stages and buffers, ``buffer_name``
+and ``stage_name``. All other integrators inherit from this one.
+
+LowStorageIntegrator
+----------------------
+
+
+The ``LowStorageIntegrator`` contains integrators in the 2S form as
+described in `Ketchson (2010)`_. These integrators are of the classic
+`Shu Osher`_ form:
+
+.. math::
+
+   u^{(0)} &= u^n \\
+   u^{(i)} &= \sum_{k=0}^{i-1} (\alpha_{i,k} u^{(k)} + \Delta t \beta_{i, k} F(u^{(k)})\\
+   u^{n+1} &= u^{(m)}
+
+where superscripts in parentheses mean subcycles in a Runge-Kutta
+integration and :math:`F` is the right-hand-side of ODE system. The
+difference between these low-storage methods and the classic SSPK
+methods is that the low-storage methods typically have sparse
+:math:`\alpha` and :math:`\beta` matrices, which are replaced by
+diagonal termes, named :math:`\gamma_0` and :math:`\gamma_1`
+respectively. 
+
+These methods can be generalized to support more general methods with
+the introduction of an additional :math:`\delta` term for first
+averaging the updated stage with previous stages. This form is also described in Section 3.2.3 of the `Athena++ paper`_.
+
+The full update then takes the form:
+
+.. math::
+
+   u^{(1)} &:= u^{(1)} + \delta_s u^{(0)} \\
+   u^{(0)} &:= \gamma_{s0} u^{(0)} + \gamma_{s1} u^{(1)} + \beta_{s,s-1} \Delta t F(u^{(0)})
+
+where here :math:`u^{(0)}` and :math:`u^{(1)}` are the two storage
+buffers required to compute the update for a given Runge-Kutta stage
+:math:`s`.
+
+.. _Ketchson (2010): https://doi.org/10.1016/j.jcp.2009.11.006
+
+.. _Shu Osher: https://doi.org/10.1016/0021-9991(88)90177-5
+
+.. _Athena++ paper: https://doi.org/10.3847/1538-4365/ab929b
+
+The ``LowStorageIntegrator`` contains arrays for ``delta``, ``beta``,
+``gam0``, and ``gam1``. Available integration methods are:
+
+* ``RK1``, which is simply forward Euler.
+
+* ``RK2``, which is Heun's method.
+
+* ``VL2``, 2nd-order Van Leer predictor-corrector from Stone and
+  Gardiner 2009. Requires donor-cell reconstruction for the predictor
+  stage.
+
+* ``RK3``, a strong stability preserving variant.
+
+* ``RK4``, a strong stability preserving variant.
+
+ButcherIntegrator
+---------------------
+
+The ``ButcherIntegrator`` provides a classic Butcher Tableau with
+arrays :math:`a` to compute the stages, :math:`c` to compute the time
+offsets, and :math:`b` to compute the final update for a time
+step. Available integration methods are:
+
+* ``RK1``, simple forward Euler.
+
+* ``RK2``, Heun's method.
+
+* ``RK4``, The classic 4th-order method.
+
+* ``RK10``, A recent version with fewer stages than Fehlberg's classic RK8(9), computed by Faegin and tabulated `here <https://sce.uhcl.edu/rungekutta/>`__.

example/particle_leapfrog/particle_leapfrog.hpp

@@ -39,7 +39,7 @@ class ParticleDriver : public EvolutionDriver {
   TaskListStatus Step();
 
  private:
-  StagedIntegrator integrator;
+  LowStorageIntegrator integrator;
 };
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin);

example/particles/particles.hpp

@@ -60,7 +60,7 @@ class ParticleDriver : public EvolutionDriver {
   TaskListStatus Step();
 
  private:
-  StagedIntegrator integrator;
+  LowStorageIntegrator integrator;
 };
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin);

src/CMakeLists.txt

@@ -202,6 +202,11 @@ add_library(parthenon
   tasks/task_list.hpp
   tasks/task_types.hpp
 
+  time_integration/butcher_integrator.cpp
+  time_integration/low_storage_integrator.cpp
+  time_integration/staged_integrator.cpp
+  time_integration/staged_integrator.hpp
+
   utils/alias_method.cpp
   utils/alias_method.hpp
   utils/buffer_utils.cpp

src/driver/multistage.cpp

@@ -16,102 +16,11 @@
 
 namespace parthenon {
 
-StagedIntegrator::StagedIntegrator(ParameterInput *pin) {
-  std::string integrator_name =
-      pin->GetOrAddString("parthenon/time", "integrator", "rk2");
-
-  if (integrator_name == "rk1") {
-    nstages = 1;
-    beta.resize(nstages);
-    gam0.resize(nstages);
-    gam1.resize(nstages);
-
-    beta[0] = 1.0;
-    gam0[0] = 0.0;
-    gam1[0] = 1.0;
-  } else if (integrator_name == "rk2") {
-    nstages = 2;
-    beta.resize(nstages);
-    gam0.resize(nstages);
-    gam1.resize(nstages);
-
-    beta[0] = 1.0;
-    gam0[0] = 0.0;
-    gam1[0] = 1.0;
-
-    beta[1] = 0.5;
-    gam0[1] = 0.5;
-    gam1[1] = 0.5;
-  } else if (integrator_name == "vl2") {
-    nstages = 2;
-    beta.resize(nstages);
-    gam0.resize(nstages);
-    gam1.resize(nstages);
-
-    beta[0] = 0.5;
-    gam0[0] = 0.0;
-    gam1[0] = 1.0;
-
-    beta[1] = 1.0;
-    gam0[1] = 0.0;
-    gam1[1] = 1.0;
-  } else if (integrator_name == "rk3") {
-    nstages = 3;
-    beta.resize(nstages);
-    gam0.resize(nstages);
-    gam1.resize(nstages);
-
-    beta[0] = 1.0;
-    gam0[0] = 0.0;
-    gam1[0] = 1.0;
-
-    beta[1] = 0.25;
-    gam0[1] = 0.25;
-    gam1[1] = 0.75;
-
-    beta[2] = 2.0 / 3.0;
-    gam0[2] = 2.0 / 3.0;
-    gam1[2] = 1.0 / 3.0;
-  } else {
-    throw std::invalid_argument("Invalid selection for the time integrator: " +
-                                integrator_name);
-  }
-  stage_name.resize(nstages + 1);
-  stage_name[0] = "base";
-  for (int i = 1; i < nstages; i++) {
-    stage_name[i] = std::to_string(i);
-  }
-  stage_name[nstages] = stage_name[0];
-}
-
-TaskListStatus MultiStageDriver::Step() {
-  Kokkos::Profiling::pushRegion("MultiStage_Step");
-  using DriverUtils::ConstructAndExecuteTaskLists;
-  TaskListStatus status;
-  integrator->dt = tm.dt;
-  for (int stage = 1; stage <= integrator->nstages; stage++) {
-    // Clear any initialization info. We should be relying
-    // on only the immediately preceding stage to contain
-    // reasonable data
-    pmesh->SetAllVariablesToInitialized();
-    status = ConstructAndExecuteTaskLists<>(this, stage);
-    if (status != TaskListStatus::complete) break;
-  }
-  Kokkos::Profiling::popRegion(); // MultiStage_Step
-  return status;
-}
-
-TaskListStatus MultiStageBlockTaskDriver::Step() {
-  Kokkos::Profiling::pushRegion("MultiStageBlockTask_Step");
-  using DriverUtils::ConstructAndExecuteBlockTasks;
-  TaskListStatus status;
-  integrator->dt = tm.dt;
-  for (int stage = 1; stage <= integrator->nstages; stage++) {
-    status = ConstructAndExecuteBlockTasks<>(this, stage);
-    if (status != TaskListStatus::complete) break;
-  }
-  Kokkos::Profiling::popRegion(); // MultiStageBlockTask_Step
-  return status;
-}
+template class MultiStageDriverGeneric<StagedIntegrator>;
+template class MultiStageBlockTaskDriverGeneric<StagedIntegrator>;
+template class MultiStageDriverGeneric<LowStorageIntegrator>;
+template class MultiStageBlockTaskDriverGeneric<LowStorageIntegrator>;
+template class MultiStageDriverGeneric<ButcherIntegrator>;
+template class MultiStageBlockTaskDriverGeneric<ButcherIntegrator>;
 
 } // namespace parthenon

src/driver/multistage.hpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2022. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -23,42 +23,67 @@
 #include "mesh/mesh.hpp"
 #include "parameter_input.hpp"
 #include "tasks/task_list.hpp"
+#include "time_integration/staged_integrator.hpp"
 
 namespace parthenon {
 
-struct StagedIntegrator {
-  StagedIntegrator() = default;
-  explicit StagedIntegrator(ParameterInput *pin);
-  int nstages;
-  std::vector<Real> beta;
-  std::vector<Real> gam0;
-  std::vector<Real> gam1;
-  std::vector<std::string> stage_name;
-  Real dt;
-};
-
-class MultiStageDriver : public EvolutionDriver {
+template <typename Integrator = LowStorageIntegrator>
+class MultiStageDriverGeneric : public EvolutionDriver {
  public:
-  MultiStageDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm)
-      : EvolutionDriver(pin, app_in, pm),
-        integrator(std::make_unique<StagedIntegrator>(pin)) {}
+  MultiStageDriverGeneric(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm)
+      : EvolutionDriver(pin, app_in, pm), integrator(std::make_unique<Integrator>(pin)) {}
   // An application driver that derives from this class must define this
   // function, which defines the application specific list of tasks and
   // the dependencies that must be executed.
   virtual TaskCollection MakeTaskCollection(BlockList_t &blocks, int stage) = 0;
-  virtual TaskListStatus Step();
+  virtual TaskListStatus Step() {
+    Kokkos::Profiling::pushRegion("MultiStage_Step");
+    using DriverUtils::ConstructAndExecuteTaskLists;
+    TaskListStatus status;
+    integrator->dt = tm.dt;
+    for (int stage = 1; stage <= integrator->nstages; stage++) {
+      // Clear any initialization info. We should be relying
+      // on only the immediately preceding stage to contain
+      // reasonable data
+      pmesh->SetAllVariablesToInitialized();
+      status = ConstructAndExecuteTaskLists<>(this, stage);
+      if (status != TaskListStatus::complete) break;
+    }
+    Kokkos::Profiling::popRegion(); // MultiStage_Step
+    return status;
+  }
 
  protected:
-  std::unique_ptr<StagedIntegrator> integrator;
+  std::unique_ptr<Integrator> integrator;
 };
+using MultiStageDriver = MultiStageDriverGeneric<LowStorageIntegrator>;
 
-class MultiStageBlockTaskDriver : public MultiStageDriver {
+template <typename Integrator = LowStorageIntegrator>
+class MultiStageBlockTaskDriverGeneric : public MultiStageDriverGeneric<Integrator> {
  public:
-  MultiStageBlockTaskDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm)
-      : MultiStageDriver(pin, app_in, pm) {}
+  MultiStageBlockTaskDriverGeneric(ParameterInput *pin, ApplicationInput *app_in,
+                                   Mesh *pm)
+      : MultiStageDriverGeneric<Integrator>(pin, app_in, pm) {}
   virtual TaskList MakeTaskList(MeshBlock *pmb, int stage) = 0;
-  virtual TaskListStatus Step();
+  virtual TaskListStatus Step() {
+    Kokkos::Profiling::pushRegion("MultiStageBlockTask_Step");
+    using DriverUtils::ConstructAndExecuteBlockTasks;
+    TaskListStatus status;
+    Integrator *integrator = (this->integrator).get();
+    SimTime tm = this->tm;
+    integrator->dt = tm.dt;
+    for (int stage = 1; stage <= integrator->nstages; stage++) {
+      status = ConstructAndExecuteBlockTasks<>(this, stage);
+      if (status != TaskListStatus::complete) break;
+    }
+    Kokkos::Profiling::popRegion(); // MultiStageBlockTask_Step
+    return status;
+  }
+
+ protected:
+  std::unique_ptr<Integrator> integrator;
 };
+using MultiStageBlockTaskDriver = MultiStageBlockTaskDriverGeneric<LowStorageIntegrator>;
 
 } // namespace parthenon