diff --git a/.gitignore b/.gitignore index 03b6f9fbbd1d..135d76be3f15 100644 --- a/.gitignore +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 7fea3baabd9e..5dfada1c7102 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/doc/sphinx/src/driver.rst b/doc/sphinx/src/driver.rst index 24b42fa09869..a7c400bc9086 100644 --- a/doc/sphinx/src/driver.rst +++ b/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 `__. 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`` object which includes members for +it in two important ways. First, it is templated on an ``Integrator`` type and +holds a ``std::unique_ptr`` 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 diff --git a/doc/sphinx/src/integrators.rst b/doc/sphinx/src/integrators.rst new file mode 100644 index 000000000000..8e035832f473 --- /dev/null +++ b/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 `__. diff --git a/example/particle_leapfrog/particle_leapfrog.hpp b/example/particle_leapfrog/particle_leapfrog.hpp index 5e2f4b478ff3..dbe910afd721 100644 --- a/example/particle_leapfrog/particle_leapfrog.hpp +++ b/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); diff --git a/example/particles/particles.hpp b/example/particles/particles.hpp index fb8c456ed30a..ca0481bb8d67 100644 --- a/example/particles/particles.hpp +++ b/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); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c5ef598442f3..fb46cad9637c 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/driver/multistage.cpp b/src/driver/multistage.cpp index 181090b0c2b6..4b81bd766a5d 100644 --- a/src/driver/multistage.cpp +++ b/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; +template class MultiStageBlockTaskDriverGeneric; +template class MultiStageDriverGeneric; +template class MultiStageBlockTaskDriverGeneric; +template class MultiStageDriverGeneric; +template class MultiStageBlockTaskDriverGeneric; } // namespace parthenon diff --git a/src/driver/multistage.hpp b/src/driver/multistage.hpp index 5935278a5202..07748ae3564c 100644 --- a/src/driver/multistage.hpp +++ b/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 beta; - std::vector gam0; - std::vector gam1; - std::vector stage_name; - Real dt; -}; - -class MultiStageDriver : public EvolutionDriver { +template +class MultiStageDriverGeneric : public EvolutionDriver { public: - MultiStageDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm) - : EvolutionDriver(pin, app_in, pm), - integrator(std::make_unique(pin)) {} + MultiStageDriverGeneric(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm) + : EvolutionDriver(pin, app_in, pm), integrator(std::make_unique(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 integrator; + std::unique_ptr integrator; }; +using MultiStageDriver = MultiStageDriverGeneric; -class MultiStageBlockTaskDriver : public MultiStageDriver { +template +class MultiStageBlockTaskDriverGeneric : public MultiStageDriverGeneric { public: - MultiStageBlockTaskDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm) - : MultiStageDriver(pin, app_in, pm) {} + MultiStageBlockTaskDriverGeneric(ParameterInput *pin, ApplicationInput *app_in, + Mesh *pm) + : MultiStageDriverGeneric(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; }; +using MultiStageBlockTaskDriver = MultiStageBlockTaskDriverGeneric; } // namespace parthenon diff --git a/src/interface/update.hpp b/src/interface/update.hpp index fec29a54c512..6a16f8b5b3ba 100644 --- a/src/interface/update.hpp +++ b/src/interface/update.hpp @@ -28,6 +28,7 @@ #include "interface/params.hpp" #include "interface/sparse_pack.hpp" #include "interface/state_descriptor.hpp" +#include "time_integration/staged_integrator.hpp" #include "kokkos_abstraction.hpp" #include "mesh/domain.hpp" @@ -87,6 +88,27 @@ TaskStatus WeightedSumData(const F &flags, T *in1, T *in2, const Real w1, const return TaskStatus::complete; } +template +TaskStatus CopyData(const F &flags, T *in, T *out) { + return WeightedSumData(flags, in, in, 1, 0, out); +} + +template +TaskStatus SetDataToConstant(const F &flags, T *data, const Real val) { + Kokkos::Profiling::pushRegion("Task_SetDataToConstant"); + const auto &x = data->PackVariables(flags); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SetDataToConstant", DevExecSpace(), 0, x.GetDim(5) - 1, 0, + x.GetDim(4) - 1, 0, x.GetDim(3) - 1, 0, x.GetDim(2) - 1, 0, x.GetDim(1) - 1, + KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + if (x.IsAllocated(b, l)) { + x(b, l, k, j, i) = val; + } + }); + Kokkos::Profiling::popRegion(); // Task_SetDataToConstant + return TaskStatus::complete; +} + template TaskStatus SumData(const F &flags, T *in1, T *in2, T *out) { return WeightedSumData(flags, in1, in2, 1.0, 1.0, out); @@ -114,6 +136,139 @@ TaskStatus AverageIndependentData(T *c1, T *c2, const Real wgt1) { (1.0 - wgt1), c1); } +// See equation 14 in Ketcheson, Jcomp 229 (2010) 1763-1773 +// In Parthenon language, s0 is the variable we are updating +// and rhs should be computed with respect to s0. +// if update_s1, s1 should be set at the beginning of the cycle to 0 +// otherwise, s1 should be set at the beginning of the RK update to be +// a copy of base. in the final stage, base for the next cycle should +// be set to s0. +template +TaskStatus Update2S(const F &flags, T *s0_data, T *s1_data, T *rhs_data, + const LowStorageIntegrator *pint, Real dt, int stage, + bool update_s1) { + Kokkos::Profiling::pushRegion("Task_2S_Update"); + const auto &s0 = s0_data->PackVariables(flags); + const auto &s1 = s1_data->PackVariables(flags); + const auto &rhs = rhs_data->PackVariables(flags); + + const IndexDomain interior = IndexDomain::interior; + const IndexRange ib = s0_data->GetBoundsI(interior); + const IndexRange jb = s0_data->GetBoundsJ(interior); + const IndexRange kb = s0_data->GetBoundsK(interior); + + Real delta = pint->delta[stage - 1]; + Real beta = pint->beta[stage - 1]; + Real gam0 = pint->gam0[stage - 1]; + Real gam1 = pint->gam1[stage - 1]; + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "2S_Update", DevExecSpace(), 0, s0.GetDim(5) - 1, 0, + s0.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + if (s0.IsAllocated(b, l) && s1.IsAllocated(b, l) && rhs.IsAllocated(b, l)) { + if (update_s1) { + s1(b, l, k, j, i) = s1(b, l, k, j, i) + delta * s0(b, l, k, j, i); + } + s0(b, l, k, j, i) = gam0 * s0(b, l, k, j, i) + gam1 * s1(b, l, k, j, i) + + beta * dt * rhs(b, l, k, j, i); + } + }); + Kokkos::Profiling::popRegion(); // Task_2S_Update + return TaskStatus::complete; +} +template +TaskStatus Update2SIndependent(T *s0_data, T *s1_data, T *rhs_data, + const LowStorageIntegrator *pint, Real dt, int stage, + bool update_s1) { + return Update2S(std::vector({Metadata::Independent}), s0_data, s1_data, + rhs_data, pint, dt, stage, update_s1); +} + +// For integration with Butcher tableaus +// returns base + dt * sum_{j=0}^{k-1} a_{kj} S_j +// for stages S_j +// This can then be used to compute right-hand sides. +template +TaskStatus SumButcher(const F &flags, std::shared_ptr base_data, + std::vector> stage_data, + std::shared_ptr out_data, const ButcherIntegrator *pint, Real dt, + int stage) { + Kokkos::Profiling::pushRegion("Task_Butcher_Sum"); + const auto &out = out_data->PackVariables(flags); + const auto &in = base_data->PackVariables(flags); + const IndexDomain interior = IndexDomain::interior; + const IndexRange ib = out_data->GetBoundsI(interior); + const IndexRange jb = out_data->GetBoundsJ(interior); + const IndexRange kb = out_data->GetBoundsK(interior); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ButcherSumInit", DevExecSpace(), 0, out.GetDim(5) - 1, 0, + out.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + if (out.IsAllocated(b, l) && in.IsAllocated(b, l)) { + out(b, l, k, j, i) = in(b, l, k, j, i); + } + }); + for (int prev = 0; prev < stage; ++prev) { + Real a = pint->a[stage - 1][prev]; + const auto &in = stage_data[stage]->PackVariables(flags); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ButcherSum", DevExecSpace(), 0, out.GetDim(5) - 1, 0, + out.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + if (out.IsAllocated(b, l) && in.IsAllocated(b, l)) { + out(b, l, k, j, i) += dt * a * in(b, l, k, j, i); + } + }); + } + Kokkos::Profiling::popRegion(); // Task_Butcher_Sum + return TaskStatus::complete; +} +template +TaskStatus SumButcherIndependent(std::shared_ptr base_data, + std::vector> stage_data, + std::shared_ptr out_data, + const ButcherIntegrator *pint, Real dt, int stage) { + return SumButcher(std::vector({Metadata::Independent}), base_data, + stage_data, out_data, pint, dt, stage); +} + +// The actual butcher update at the final stage of a cycle +template +TaskStatus UpdateButcher(const F &flags, std::vector> stage_data, + std::shared_ptr out_data, const ButcherIntegrator *pint, + Real dt) { + Kokkos::Profiling::pushRegion("Task_Butcher_Update"); + + const auto &out = out_data->PackVariables(flags); + const IndexDomain interior = IndexDomain::interior; + const IndexRange ib = out_data->GetBoundsI(interior); + const IndexRange jb = out_data->GetBoundsJ(interior); + const IndexRange kb = out_data->GetBoundsK(interior); + + const int nstages = pint->nstages; + for (int stage = 0; stage < nstages; ++stage) { + const Real butcher_b = pint->b[stage]; + const auto &in = stage_data[stage]->PackVariables(flags); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ButcherUpdate", DevExecSpace(), 0, out.GetDim(5) - 1, 0, + out.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int l, const int k, const int j, const int i) { + if (out.IsAllocated(b, l) && in.IsAllocated(b, l)) { + out(b, l, k, j, i) += dt * b * in(b, l, k, j, i); + } + }); + } + Kokkos::Profiling::popRegion(); // Task_Butcher_Update + return TaskStatus::complete; +} +template +TaskStatus UpdateButcherIndependent(std::vector> stage_data, + std::shared_ptr out_data, + const ButcherIntegrator *pint, Real dt) { + return UpdateButcherIndependent(std::vector({Metadata::Independent}), + stage_data, out_data, pint, dt); +} + template TaskStatus EstimateTimestep(T *rc) { Kokkos::Profiling::pushRegion("Task_EstimateTimestep"); diff --git a/src/time_integration/butcher_integrator.cpp b/src/time_integration/butcher_integrator.cpp new file mode 100644 index 000000000000..d90de479ec68 --- /dev/null +++ b/src/time_integration/butcher_integrator.cpp @@ -0,0 +1,231 @@ +//======================================================================================== +// (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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include + +#include "basic_types.hpp" +#include "parameter_input.hpp" +#include "staged_integrator.hpp" + +namespace parthenon { + +/* + *classic Butcher Tableau integrators + * + * Notation: + * c | a + * --+-- + * | b + * + * Note that a good resource is NASA-TR-R-287 + * https://ntrs.nasa.gov/citations/19680027281 + * Note that in the report, Fehlberg used the following notation, + * as translated into standard Butcher tableau notation: + * beta_{k lambda} = a_{ij} + * alpha_k = c_k + * c_k = b_k + */ +ButcherIntegrator::ButcherIntegrator(ParameterInput *pin) + : StagedIntegrator(pin->GetOrAddString("parthenon/time", "integrator", "rk2")) { + if (name_ == "rk1") { + nstages = nbuffers = 1; + Resize_(nstages); + + a[0][0] = 0; + b[0] = 1; + c[0] = 0; + } else if (name_ == "rk2") { + // Heun's method. Should match minimal storage solution + nstages = nbuffers = 2; + Resize_(nstages); + + a[0] = {0, 0}; + a[1] = {1, 0}; + b = {0, 1}; + c = {0, 1. / 3., 2. / 3.}; + } else if (name_ == "rk4") { + // Classic RK4 because why not + nstages = nbuffers = 4; + Resize_(nstages); + + /* clang-format off */ + a[0] = {0, 0, 0, 0}; + a[1] = {0.5, 0, 0, 0}; + a[2] = {0, 0.5, 0, 0}; + a[3] = {0, 0, 1, 0}; + /* clang-format on */ + b = {1. / 6., 1. / 3., 1. / 3., 1. / 6.}; + c = {0, 0.5, 0.5, 1}; + } else if (name_ == "rk10") { + // Feagin's family of high-order embedded methods as introduced in + // Feagin, Neural, Parallel, and Scientific Computations 20 (2012) + // 437-458 + // Note that Feagin uses the following notation, which we have + // translated: + // + // a | beta + // --+----- + // | c + // + // Feagin's 10th-order method has the same number of stages + // as Fehlberg's classic RK8(9) adaptive method + // + // Feagin's coefficients are vailable on his website: + // https://sce.uhcl.edu/rungekutta/ + + nstages = nbuffers = 17; + Resize_(nstages); + + // computed up to 60 digits + c[0] = 0.000000000000000000000000000000000000000000000000000000000000; + c[1] = 0.100000000000000000000000000000000000000000000000000000000000; + c[2] = 0.539357840802981787532485197881302436857273449701009015505500; + c[3] = 0.809036761204472681298727796821953655285910174551513523258250; + c[4] = 0.309036761204472681298727796821953655285910174551513523258250; + c[5] = 0.981074190219795268254879548310562080489056746118724882027805; + c[6] = 0.833333333333333333333333333333333333333333333333333333333333; + c[7] = 0.354017365856802376329264185948796742115824053807373968324184; + c[8] = 0.882527661964732346425501486979669075182867844268052119663791; + c[9] = 0.642615758240322548157075497020439535959501736363212695909875; + c[10] = 0.357384241759677451842924502979560464040498263636787304090125; + c[11] = 0.117472338035267653574498513020330924817132155731947880336209; + c[12] = 0.833333333333333333333333333333333333333333333333333333333333; + c[13] = 0.309036761204472681298727796821953655285910174551513523258250; + c[14] = 0.539357840802981787532485197881302436857273449701009015505500; + c[15] = 0.100000000000000000000000000000000000000000000000000000000000; + c[16] = 1.00000000000000000000000000000000000000000000000000000000000; + + b[0] = 0.0333333333333333333333333333333333333333333333333333333333333; + b[1] = 0.0250000000000000000000000000000000000000000000000000000000000; + b[2] = 0.0333333333333333333333333333333333333333333333333333333333333; + b[3] = 0.000000000000000000000000000000000000000000000000000000000000; + b[4] = 0.0500000000000000000000000000000000000000000000000000000000000; + b[5] = 0.000000000000000000000000000000000000000000000000000000000000; + b[6] = 0.0400000000000000000000000000000000000000000000000000000000000; + b[7] = 0.000000000000000000000000000000000000000000000000000000000000; + b[8] = 0.189237478148923490158306404106012326238162346948625830327194; + b[9] = 0.277429188517743176508360262560654340428504319718040836339472; + b[10] = 0.277429188517743176508360262560654340428504319718040836339472; + b[11] = 0.189237478148923490158306404106012326238162346948625830327194; + b[12] = -0.0400000000000000000000000000000000000000000000000000000000000; + b[13] = -0.0500000000000000000000000000000000000000000000000000000000000; + b[14] = -0.0333333333333333333333333333333333333333333333333333333333333; + b[15] = -0.0250000000000000000000000000000000000000000000000000000000000; + b[16] = 0.0333333333333333333333333333333333333333333333333333333333333; + + // Set matrix to zero first, to only set non-zero coeffs + for (int i = 0; i < nstages; ++i) { + for (int j = 0; j < nstages; ++j) { + a[i][j] = 0; + } + } + a[1][0] = 0.100000000000000000000000000000000000000000000000000000000000; + a[2][0] = -0.915176561375291440520015019275342154318951387664369720564660; + a[2][1] = 1.45453440217827322805250021715664459117622483736537873607016; + a[3][0] = 0.202259190301118170324681949205488413821477543637878380814562; + a[3][2] = 0.606777570903354510974045847616465241464432630913635142443687; + a[4][0] = 0.184024714708643575149100693471120664216774047979591417844635; + a[4][2] = 0.197966831227192369068141770510388793370637287463360401555746; + a[4][3] = -0.0729547847313632629185146671595558023015011608914382961421311; + a[5][0] = 0.0879007340206681337319777094132125475918886824944548534041378; + a[5][3] = 0.410459702520260645318174895920453426088035325902848695210406; + a[5][4] = 0.482713753678866489204726942976896106809132737721421333413261; + a[6][0] = 0.0859700504902460302188480225945808401411132615636600222593880; + a[6][3] = 0.330885963040722183948884057658753173648240154838402033448632; + a[6][4] = 0.489662957309450192844507011135898201178015478433790097210790; + a[6][5] = -0.0731856375070850736789057580558988816340355615025188195854775; + a[7][0] = 0.120930449125333720660378854927668953958938996999703678812621; + a[7][4] = 0.260124675758295622809007617838335174368108756484693361887839; + a[7][5] = 0.0325402621549091330158899334391231259332716675992700000776101; + a[7][6] = -0.0595780211817361001560122202563305121444953672762930724538856; + a[8][0] = 0.110854379580391483508936171010218441909425780168656559807038; + a[8][5] = -0.0605761488255005587620924953655516875526344415354339234619466; + a[8][6] = 0.321763705601778390100898799049878904081404368603077129251110; + a[8][7] = 0.510485725608063031577759012285123416744672137031752354067590; + a[9][0] = 0.112054414752879004829715002761802363003717611158172229329393; + a[9][5] = -0.144942775902865915672349828340980777181668499748506838876185; + a[9][6] = -0.333269719096256706589705211415746871709467423992115497968724; + a[9][7] = 0.499269229556880061353316843969978567860276816592673201240332; + a[9][8] = 0.509504608929686104236098690045386253986643232352989602185060; + a[10][0] = 0.113976783964185986138004186736901163890724752541486831640341; + a[10][5] = -0.0768813364203356938586214289120895270821349023390922987406384; + a[10][6] = 0.239527360324390649107711455271882373019741311201004119339563; + a[10][7] = 0.397774662368094639047830462488952104564716416343454639902613; + a[10][8] = 0.0107558956873607455550609147441477450257136782823280838547024; + a[10][9] = -0.327769124164018874147061087350233395378262992392394071906457; + a[11][0] = 0.0798314528280196046351426864486400322758737630423413945356284; + a[11][5] = -0.0520329686800603076514949887612959068721311443881683526937298; + a[11][6] = -0.0576954146168548881732784355283433509066159287152968723021864; + a[11][7] = 0.194781915712104164976306262147382871156142921354409364738090; + a[11][8] = 0.145384923188325069727524825977071194859203467568236523866582; + a[11][9] = -0.0782942710351670777553986729725692447252077047239160551335016; + a[11][10] = -0.114503299361098912184303164290554670970133218405658122674674; + a[12][0] = 0.985115610164857280120041500306517278413646677314195559520529; + a[12][3] = 0.330885963040722183948884057658753173648240154838402033448632; + a[12][4] = 0.489662957309450192844507011135898201178015478433790097210790; + a[12][5] = -1.37896486574843567582112720930751902353904327148559471526397; + a[12][6] = -0.861164195027635666673916999665534573351026060987427093314412; + a[12][7] = 5.78428813637537220022999785486578436006872789689499172601856; + a[12][8] = 3.28807761985103566890460615937314805477268252903342356581925; + a[12][9] = -2.38633905093136384013422325215527866148401465975954104585807; + a[12][10] = -3.25479342483643918654589367587788726747711504674780680269911; + a[12][11] = -2.16343541686422982353954211300054820889678036420109999154887; + a[13][0] = 0.895080295771632891049613132336585138148156279241561345991710; + a[13][2] = 0.197966831227192369068141770510388793370637287463360401555746; + a[13][3] = -0.0729547847313632629185146671595558023015011608914382961421311; + a[13][5] = -0.851236239662007619739049371445966793289359722875702227166105; + a[13][6] = 0.398320112318533301719718614174373643336480918103773904231856; + a[13][7] = 3.63937263181035606029412920047090044132027387893977804176229; + a[13][8] = 1.54822877039830322365301663075174564919981736348973496313065; + a[13][9] = -2.12221714704053716026062427460427261025318461146260124401561; + a[13][10] = -1.58350398545326172713384349625753212757269188934434237975291; + a[13][11] = -1.71561608285936264922031819751349098912615880827551992973034; + a[13][12] = -0.0244036405750127452135415444412216875465593598370910566069132; + a[14][0] = -0.915176561375291440520015019275342154318951387664369720564660; + a[14][1] = 1.45453440217827322805250021715664459117622483736537873607016; + a[14][4] = -0.777333643644968233538931228575302137803351053629547286334469; + a[14][6] = -0.0910895662155176069593203555807484200111889091770101799647985; + a[14][12] = 0.0910895662155176069593203555807484200111889091770101799647985; + a[14][13] = 0.777333643644968233538931228575302137803351053629547286334469; + a[15][0] = 0.100000000000000000000000000000000000000000000000000000000000; + a[15][2] = -0.157178665799771163367058998273128921867183754126709419409654; + a[15][14] = 0.157178665799771163367058998273128921867183754126709419409654; + a[16][0] = 0.181781300700095283888472062582262379650443831463199521664945; + a[16][1] = 0.675000000000000000000000000000000000000000000000000000000000; + a[16][2] = 0.342758159847189839942220553413850871742338734703958919937260; + a[16][4] = 0.259111214548322744512977076191767379267783684543182428778156; + a[16][5] = -0.358278966717952089048961276721979397739750634673268802484271; + a[16][6] = -1.04594895940883306095050068756409905131588123172378489286080; + a[16][7] = 0.930327845415626983292300564432428777137601651182965794680397; + a[16][8] = 1.77950959431708102446142106794824453926275743243327790536000; + a[16][9] = 0.100000000000000000000000000000000000000000000000000000000000; + a[16][10] = -0.282547569539044081612477785222287276408489375976211189952877; + a[16][11] = -0.159327350119972549169261984373485859278031542127551931461821; + a[16][12] = -0.145515894647001510860991961081084111308650130578626404945571; + a[16][13] = -0.259111214548322744512977076191767379267783684543182428778156; + a[16][14] = -0.342758159847189839942220553413850871742338734703958919937260; + a[16][15] = -0.675000000000000000000000000000000000000000000000000000000000; + } +} + +void ButcherIntegrator::Resize_(int nstages) { + a.resize(nstages); + for (int i = 0; i < a.size(); ++i) { + a[i].resize(nstages); + } + b.resize(nstages); + c.resize(nstages); +} + +} // namespace parthenon diff --git a/src/time_integration/low_storage_integrator.cpp b/src/time_integration/low_storage_integrator.cpp new file mode 100644 index 000000000000..ab874d98e050 --- /dev/null +++ b/src/time_integration/low_storage_integrator.cpp @@ -0,0 +1,158 @@ +//======================================================================================== +// (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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include + +#include "basic_types.hpp" +#include "parameter_input.hpp" +#include "staged_integrator.hpp" + +namespace parthenon { + +/* + * These integrators are of the 2S form as described in + * Ketcheson, Jcomp 229 (2010) 1763-1773 + * See Equation 14. + * + * These integrators are of the classic Shu Osher form + * 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) + * See Shu and Osher, JComp 77 (1988) 439-471 + * + * The difference between these low-storage methods and classic SSPK + * methods in Shu Osher form is that low-storage methods typically have + * sparse alpha and beta matrices, meaning fewer past stages are + * needed. Here the alpha and beta matrices are replaced by their + * diagonal terms, named gamma0 and gamma1. + * + * They can be generalized to support more general methods with the + * introduction of a delta term for a first averaging. + * The form is also described in Section 3.2.3 of the Athena++ paper: + * Stone et al., ApJS (2020) 249:4 + * See equations 11 through 15. + */ +LowStorageIntegrator::LowStorageIntegrator(ParameterInput *pin) + : StagedIntegrator(pin->GetOrAddString("parthenon/time", "integrator", "rk2")) { + if (name_ == "rk1") { + nstages = 1; + nbuffers = 1; + delta.resize(nstages); + beta.resize(nstages); + gam0.resize(nstages); + gam1.resize(nstages); + + delta[0] = 1.0; + beta[0] = 1.0; + gam0[0] = 0.0; + gam1[0] = 1.0; + } else if (name_ == "rk2") { + nstages = 2; + nbuffers = 2; + delta.resize(nstages); + beta.resize(nstages); + gam0.resize(nstages); + gam1.resize(nstages); + + delta[0] = 1.0; + beta[0] = 1.0; + gam0[0] = 0.0; + gam1[0] = 1.0; + + delta[1] = 0.0; + beta[1] = 0.5; + gam0[1] = 0.5; + gam1[1] = 0.5; + } else if (name_ == "vl2") { + nstages = 2; + nbuffers = 2; + delta.resize(nstages); + beta.resize(nstages); + gam0.resize(nstages); + gam1.resize(nstages); + + delta[0] = 1.0; + beta[0] = 0.5; + gam0[0] = 0.0; + gam1[0] = 1.0; + + delta[1] = 0.0; + beta[1] = 1.0; + gam0[1] = 0.0; + gam1[1] = 1.0; + } else if (name_ == "rk3") { + nstages = 3; + nbuffers = 2; + delta.resize(nstages); + beta.resize(nstages); + gam0.resize(nstages); + gam1.resize(nstages); + + delta[0] = 1.0; + beta[0] = 1.0; + gam0[0] = 0.0; + gam1[0] = 1.0; + + delta[1] = 0.0; + beta[1] = 0.25; + gam0[1] = 0.25; + gam1[1] = 0.75; + + delta[2] = 0.0; + beta[2] = 2.0 / 3.0; + gam0[2] = 2.0 / 3.0; + gam1[2] = 1.0 / 3.0; + } else if (name_ == "rk4") { + // Classic 5-stage SSPRK(5)4 in low-storage form + // ceff = 0.377 + // From Table 4 of Ketcheson, Jcomp 229 (2010) 1763-1773 + nstages = 5; + nbuffers = 2; + delta.resize(nstages); + beta.resize(nstages); + gam0.resize(nstages); + gam1.resize(nstages); + + delta[0] = 1.0; + beta[0] = 0.357534921136978; + gam0[0] = 0.0; + gam1[0] = 1.0; + + delta[1] = 0.0; + beta[1] = 2.364680399061355; + gam0[1] = -3.666545952121251; + gam1[1] = 4.666545952121251; + + delta[2] = 0.0; + beta[2] = 0.016239790859612; + gam0[2] = 0.035802535958088; + gam1[2] = 0.964197464041912; + + delta[3] = 0.0; + beta[3] = 0.498173799587251; + gam0[3] = 4.398279365655791; + gam1[3] = -3.398279365655790; + + delta[4] = 0.0; + beta[4] = 0.433334235669763; + gam0[4] = 0.770411587328417; + gam1[4] = 0.229588412671583; + } else { + throw std::invalid_argument("Invalid selection for the time integrator: " + name_); + } + MakePeriodicNames_(buffer_name, nbuffers); + MakePeriodicNames_(stage_name, nstages); +} + +} // namespace parthenon diff --git a/src/time_integration/staged_integrator.cpp b/src/time_integration/staged_integrator.cpp new file mode 100644 index 000000000000..22b96342703a --- /dev/null +++ b/src/time_integration/staged_integrator.cpp @@ -0,0 +1,32 @@ +//======================================================================================== +// (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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include + +#include "basic_types.hpp" +#include "parameter_input.hpp" +#include "staged_integrator.hpp" + +namespace parthenon { + +void StagedIntegrator::MakePeriodicNames_(std::vector &names, int n) { + names.resize(n + 1); + names[0] = "base"; + for (int i = 1; i < n; i++) { + names[i] = std::to_string(i); + } + names[n] = names[0]; +} + +} // namespace parthenon diff --git a/src/time_integration/staged_integrator.hpp b/src/time_integration/staged_integrator.hpp new file mode 100644 index 000000000000..0d22d7d85600 --- /dev/null +++ b/src/time_integration/staged_integrator.hpp @@ -0,0 +1,74 @@ +//======================================================================================== +// (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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef TIME_INTEGRATION_STAGED_INTEGRATOR_HPP_ +#define TIME_INTEGRATION_STAGED_INTEGRATOR_HPP_ + +#include +#include + +#include "basic_types.hpp" +#include "parameter_input.hpp" + +namespace parthenon { +class StagedIntegrator { + public: + StagedIntegrator() = default; + explicit StagedIntegrator(const std::string &name) : name_(name) {} + explicit StagedIntegrator(ParameterInput *pin) + : name_(pin->GetOrAddString("parthenon/time", "integrator", "rk2")) {} + + // JMM: These might be better as private with accessors, but I will + // leave them public for backwards compatibility. + Real dt; + int nstages; // number of stages/steps in integration + int nbuffers; // number of "time levels" of that need to be allocated + // Names of time levels. + std::vector buffer_name; + // Names of integration stages (for backwards compatibility) + // TODO(JMM): Remove this eventually + std::vector stage_name; + + const std::string &GetName() const { return name_; } + + protected: + std::string name_; + void MakePeriodicNames_(std::vector &names, int n); +}; + +class LowStorageIntegrator : public StagedIntegrator { + public: + LowStorageIntegrator() = default; + explicit LowStorageIntegrator(ParameterInput *pin); + std::vector delta; + std::vector beta; + std::vector gam0; + std::vector gam1; +}; + +// TODO(JMM): Should this be named ButcherTableauIntegrator? +class ButcherIntegrator : public StagedIntegrator { + public: + ButcherIntegrator() = default; + explicit ButcherIntegrator(ParameterInput *pin); + // TODO(JMM): Should I do a flat array with indexing instead? + std::vector> a; + std::vector b, c; + + protected: + void Resize_(int nstages); +}; + +} // namespace parthenon + +#endif // TIME_INTEGRATION_STAGED_INTEGRATOR_HPP_ diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt index 73f2414b93ae..8d1d05c5d30e 100644 --- a/tst/unit/CMakeLists.txt +++ b/tst/unit/CMakeLists.txt @@ -37,6 +37,7 @@ list(APPEND unit_tests_SOURCES test_error_checking.cpp test_partitioning.cpp test_state_descriptor.cpp + test_unit_integrators.cpp ) add_executable(unit_tests "${unit_tests_SOURCES}") diff --git a/tst/unit/test_unit_integrators.cpp b/tst/unit/test_unit_integrators.cpp new file mode 100644 index 000000000000..e6ab151e05c2 --- /dev/null +++ b/tst/unit/test_unit_integrators.cpp @@ -0,0 +1,250 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2022 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "basic_types.hpp" +#include "parameter_input.hpp" +#include "time_integration/staged_integrator.hpp" + +using parthenon::ButcherIntegrator; +using parthenon::LowStorageIntegrator; +using parthenon::ParameterInput; +using parthenon::Real; +using parthenon::StagedIntegrator; + +// Test our integrators by integrating the equation +// for a harmonic oscillator: +// +// d^2 u/dt^2 = - k^2 u +// +// reduced in 1d to: +// +// du/dt = v +// dv/dt = -k^2 u +// +// This is a strenuous test since any integrator that isn't +// time-reversible will accrue error every period. + +constexpr Real K = 2 * M_PI; +constexpr std::size_t NVARS = 2; +using State_t = std::array; +void GetRHS(const State_t &u, State_t &rhs) { + rhs[0] = u[1]; + rhs[1] = -K * K * u[0]; +} + +void GetTrueSolution(const Real t, State_t &u) { + u[0] = std::cos(K * t); + u[1] = -K * std::sin(K * t); +} + +void GetInitialData(State_t &u) { GetTrueSolution(0, u); } + +template +auto MakeIntegrator(const std::string &integration_strategy) { + ParameterInput in; + in.SetString("parthenon/time", "integrator", integration_strategy); + return T(&in); +} + +/* + * See Equation 14 in + * Ketcheson, Jcomp 229 (2010) 1763-1773 + */ +void Step2SStar(const LowStorageIntegrator &integrator, Real dt, State_t &u) { + const int nstages = integrator.nstages; + State_t &S0 = u; + State_t S1; + for (int v = 0; v < NVARS; ++v) { + S1[v] = u[v]; // set "last step" to prepare for next cycle + } + for (int stage = 1; stage <= nstages; stage++) { + State_t rhs; + const Real delta = integrator.delta[stage - 1]; + const Real beta = integrator.beta[stage - 1]; + const Real gam0 = integrator.gam0[stage - 1]; + const Real gam1 = integrator.gam1[stage - 1]; + GetRHS(S0, rhs); + for (int v = 0; v < NVARS; ++v) { + // S1[v] = S1[v] + delta * S0[v]; + S0[v] = gam0 * S0[v] + gam1 * S1[v] + beta * dt * rhs[v]; + } + } + for (int v = 0; v < NVARS; ++v) { + u[v] = S0[v]; + } +} + +void StepButcher(const ButcherIntegrator &integrator, Real dt, State_t &u) { + const int nstages = integrator.nstages; + std::vector K(nstages); + for (int stage = 0; stage < nstages; ++stage) { + State_t scratch; + for (int v = 0; v < NVARS; ++v) { + scratch[v] = u[v]; + } + for (int prev = 0; prev < stage; ++prev) { + for (int v = 0; v < NVARS; ++v) { + scratch[v] += dt * integrator.a[stage][prev] * K[prev][v]; + } + } + GetRHS(scratch, K[stage]); + } + for (int stage = 0; stage < nstages; ++stage) { + for (int v = 0; v < NVARS; ++v) { + u[v] += dt * integrator.b[stage] * K[stage][v]; + } + } +} + +template +void Integrate(const Integrator &integrator, const Stepper &step, const Real tf, Real dt, + State_t &u0) { + assert(tf > 0); + assert(dt > 0); + assert(dt < tf); + // stupid game to align to EXACTLY tf + int NT = std::ceil(tf / dt); + Real t = 0; + for (int i = 0; i < NT; ++i) { + step(integrator, dt, u0); + t += dt; + } + if ((t < tf) && (std::abs(tf - t) > 1e-12)) { + dt = tf - t; + step(integrator, dt, u0); + t += dt; + } +} + +TEST_CASE("Low storage integrator", "[StagedIntegrator]") { + GIVEN("A state with an initial condition") { + Real t0 = 0, tf = 1.15; // delibarately not a nice fraction of a period + State_t ufinal; + GetTrueSolution(tf, ufinal); + WHEN("We integrate with LowStorage rk1") { + constexpr Real dt = 1e-5; + auto integrator = MakeIntegrator("rk1"); + State_t u; + GetInitialData(u); + Integrate(integrator, Step2SStar, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with LowStorage rk2") { + constexpr Real dt = 1e-3; + auto integrator = MakeIntegrator("rk2"); + State_t u; + GetInitialData(u); + Integrate(integrator, Step2SStar, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with LowStorage vl2") { + constexpr Real dt = 1e-3; + auto integrator = MakeIntegrator("vl2"); + State_t u; + GetInitialData(u); + Integrate(integrator, Step2SStar, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with LowStorage rk3") { + constexpr Real dt = 1e-2; + auto integrator = MakeIntegrator("rk3"); + State_t u; + GetInitialData(u); + Integrate(integrator, Step2SStar, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with LowStorage rk4") { + // still accurate with large timestep + constexpr Real dt = 1e-2; + auto integrator = MakeIntegrator("rk4"); + State_t u; + GetInitialData(u); + Integrate(integrator, Step2SStar, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with butcher rk1") { + constexpr Real dt = 1e-5; + auto integrator = MakeIntegrator("rk1"); + State_t u; + GetInitialData(u); + Integrate(integrator, StepButcher, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with butcher rk2") { + constexpr Real dt = 1e-5; + auto integrator = MakeIntegrator("rk2"); + State_t u; + GetInitialData(u); + Integrate(integrator, StepButcher, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-2); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-2); + } + } + WHEN("We integrate with butcher rk4") { + // Still good accuracy with large timestep. + constexpr Real dt = 1e-2; + auto integrator = MakeIntegrator("rk4"); + State_t u; + GetInitialData(u); + Integrate(integrator, StepButcher, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-4); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-4); + } + } + WHEN("We integrate with butcher rk10") { + constexpr Real dt = 5e-2; // appears to be largest stable timestep + auto integrator = MakeIntegrator("rk10"); + State_t u; + GetInitialData(u); + Integrate(integrator, StepButcher, tf, dt, u); + THEN("The final state doesn't differ too much from the true solution") { + REQUIRE(std::abs(u[0] - ufinal[0]) <= 1e-9); + REQUIRE(std::abs(u[1] - ufinal[1]) <= 1e-9); + } + } + } +}