Add Kennedy and Carpenter's IMEX time steppers

docs/References.bib

@@ -1247,6 +1247,20 @@ @article{Kastaun2020uxr
   year =         2021
 }
 
+@article{Kennedy2003,
+  author = "Kennedy, Christopher A. and Carpenter, Mark H.",
+  title = "Additive {Runge-Kutta} schemes for convection-diffusion-reaction
+           equations",
+  journal = "Applied Numerical Mathematics",
+  volume = "44",
+  number = "1",
+  pages = "139--181",
+  year = "2003",
+  doi = "10.1016/S0168-9274(02)00138-1",
+  issn = "0168-9274",
+  url = "https://www.sciencedirect.com/science/article/pii/S0168927402001381",
+}
+
 @techreport{Kennedy2016,
   author = "Kennedy, Christopher A. and Carpenter, Mark H.",
   title = "Diagonally Implicit {Runge-Kutta} Methods for Ordinary Differential

src/Time/TimeSteppers/CMakeLists.txt

@@ -12,7 +12,9 @@ spectre_target_sources(
   Heun2.cpp
   ImexRungeKutta.cpp
   Rk3HesthavenSsp.cpp
+  Rk3Kennedy.cpp
   Rk3Owren.cpp
+  Rk4Kennedy.cpp
   Rk4Owren.cpp
   Rk5Owren.cpp
   Rk5Tsitouras.cpp
@@ -34,7 +36,9 @@ spectre_target_headers(
   ImexTimeStepper.hpp
   LtsTimeStepper.hpp
   Rk3HesthavenSsp.hpp
+  Rk3Kennedy.hpp
   Rk3Owren.hpp
+  Rk4Kennedy.hpp
   Rk4Owren.hpp
   Rk5Owren.hpp
   Rk5Tsitouras.hpp

src/Time/TimeSteppers/Factory.hpp

@@ -9,7 +9,9 @@
 #include "Time/TimeSteppers/DormandPrince5.hpp"
 #include "Time/TimeSteppers/Heun2.hpp"
 #include "Time/TimeSteppers/Rk3HesthavenSsp.hpp"
+#include "Time/TimeSteppers/Rk3Kennedy.hpp"
 #include "Time/TimeSteppers/Rk3Owren.hpp"
+#include "Time/TimeSteppers/Rk4Kennedy.hpp"
 #include "Time/TimeSteppers/Rk4Owren.hpp"
 #include "Time/TimeSteppers/Rk5Owren.hpp"
 #include "Time/TimeSteppers/Rk5Tsitouras.hpp"
@@ -21,12 +23,15 @@ using time_steppers =
     tmpl::list<TimeSteppers::AdamsBashforth, TimeSteppers::AdamsMoultonPc,
                TimeSteppers::ClassicalRungeKutta4, TimeSteppers::DormandPrince5,
                TimeSteppers::Heun2, TimeSteppers::Rk3HesthavenSsp,
-               TimeSteppers::Rk3Owren, TimeSteppers::Rk4Owren,
+               TimeSteppers::Rk3Kennedy, TimeSteppers::Rk3Owren,
+               TimeSteppers::Rk4Kennedy, TimeSteppers::Rk4Owren,
                TimeSteppers::Rk5Owren, TimeSteppers::Rk5Tsitouras>;
 
 /// Typelist of available LtsTimeSteppers
 using lts_time_steppers = tmpl::list<TimeSteppers::AdamsBashforth>;
 
 /// Typelist of available ImexTimeSteppers
-using imex_time_steppers = tmpl::list<TimeSteppers::Heun2>;
-}  // namespace Triggers
+using imex_time_steppers =
+    tmpl::list<TimeSteppers::Heun2, TimeSteppers::Rk3Kennedy,
+               TimeSteppers::Rk4Kennedy>;
+}  // namespace TimeSteppers

src/Time/TimeSteppers/Rk3Kennedy.cpp

@@ -0,0 +1,60 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Time/TimeSteppers/Rk3Kennedy.hpp"
+
+namespace TimeSteppers {
+
+size_t Rk3Kennedy::order() const { return 3; }
+
+size_t Rk3Kennedy::error_estimate_order() const { return 2; }
+
+double Rk3Kennedy::stable_step() const { return 1.832102281377816; }
+
+size_t Rk3Kennedy::imex_order() const { return 3; }
+
+size_t Rk3Kennedy::implicit_stage_order() const { return 2; }
+
+// The numbers in the tableaus are not exact.  Despite being written
+// as rationals, the true values are irrational numbers.
+const RungeKutta::ButcherTableau& Rk3Kennedy::butcher_tableau() const {
+  static const ButcherTableau tableau{
+      // Substep times
+      {1767732205903.0 / 2027836641118.0, 3.0 / 5.0, 1.0},
+      // Substep coefficients
+      {{1767732205903.0 / 2027836641118.0},
+       {5535828885825.0 / 10492691773637.0, 788022342437.0 / 10882634858940.0},
+       {6485989280629.0 / 16251701735622.0, -4246266847089.0 / 9704473918619.0,
+        10755448449292.0 / 10357097424841.0}},
+      // Result coefficients
+      {1471266399579.0 / 7840856788654.0, -4482444167858.0 / 7529755066697.0,
+       11266239266428.0 / 11593286722821.0, 1767732205903.0 / 4055673282236.0},
+      // Coefficients for the embedded method for generating an error measure.
+      {2756255671327.0 / 12835298489170.0, -10771552573575.0 / 22201958757719.0,
+       9247589265047.0 / 10645013368117.0, 2193209047091.0 / 5459859503100.0},
+      // Dense output coefficient polynomials
+      {{0.0, 4655552711362.0 / 22874653954995.0,
+        -215264564351.0 / 13552729205753.0},
+       {0.0, -18682724506714.0 / 9892148508045.0,
+        17870216137069.0 / 13817060693119.0},
+       {0.0, 34259539580243.0 / 13192909600954.0,
+        -28141676662227.0 / 17317692491321.0},
+       {0.0, 584795268549.0 / 6622622206610.0,
+        2508943948391.0 / 7218656332882.0}}};
+  return tableau;
+}
+
+const ImexRungeKutta::ImplicitButcherTableau&
+Rk3Kennedy::implicit_butcher_tableau() const {
+  static const ImplicitButcherTableau tableau{
+      {{1767732205903.0 / 4055673282236.0, 1767732205903.0 / 4055673282236.0},
+       {2746238789719.0 / 10658868560708.0, -640167445237.0 / 6845629431997.0,
+        1767732205903.0 / 4055673282236.0},
+       {1471266399579.0 / 7840856788654.0, -4482444167858.0 / 7529755066697.0,
+        11266239266428.0 / 11593286722821.0,
+        1767732205903.0 / 4055673282236.0}}};
+  return tableau;
+}
+}  // namespace TimeSteppers
+
+PUP::able::PUP_ID TimeSteppers::Rk3Kennedy::my_PUP_ID = 0;  // NOLINT

src/Time/TimeSteppers/Rk3Kennedy.hpp

@@ -0,0 +1,64 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "Options/String.hpp"
+#include "Time/TimeSteppers/ImexRungeKutta.hpp"
+#include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace TimeSteppers {
+/*!
+ * \ingroup TimeSteppersGroup
+ * \brief A third-order Runge-Kutta method with IMEX support.
+ *
+ * The coefficients are given as ARK3(2)4L[2]SA in \cite Kennedy2003.
+ *
+ * The implicit part is stiffly accurate and L-stable.
+ *
+ * The CFL factor/stable step size is 1.832102281377816.
+ */
+class Rk3Kennedy : public ImexRungeKutta {
+ public:
+  using options = tmpl::list<>;
+  static constexpr Options::String help = {
+      "A 3rd-order Runge-Kutta scheme devised by Kennedy and Carpenter."};
+
+  Rk3Kennedy() = default;
+  Rk3Kennedy(const Rk3Kennedy&) = default;
+  Rk3Kennedy& operator=(const Rk3Kennedy&) = default;
+  Rk3Kennedy(Rk3Kennedy&&) = default;
+  Rk3Kennedy& operator=(Rk3Kennedy&&) = default;
+  ~Rk3Kennedy() override = default;
+
+  size_t order() const override;
+
+  size_t error_estimate_order() const override;
+
+  double stable_step() const override;
+
+  size_t imex_order() const override;
+
+  size_t implicit_stage_order() const override;
+
+  WRAPPED_PUPable_decl_template(Rk3Kennedy);  // NOLINT
+
+  explicit Rk3Kennedy(CkMigrateMessage* /*unused*/) {}
+
+  const ButcherTableau& butcher_tableau() const override;
+
+  const ImplicitButcherTableau& implicit_butcher_tableau() const override;
+};
+
+inline bool constexpr operator==(const Rk3Kennedy& /*lhs*/,
+                                 const Rk3Kennedy& /*rhs*/) {
+  return true;
+}
+
+inline bool constexpr operator!=(const Rk3Kennedy& lhs, const Rk3Kennedy& rhs) {
+  return not(lhs == rhs);
+}
+}  // namespace TimeSteppers

src/Time/TimeSteppers/Rk4Kennedy.cpp

@@ -0,0 +1,71 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Time/TimeSteppers/Rk4Kennedy.hpp"
+
+namespace TimeSteppers {
+
+size_t Rk4Kennedy::order() const { return 4; }
+
+size_t Rk4Kennedy::error_estimate_order() const { return 3; }
+
+double Rk4Kennedy::stable_step() const { return 2.1172491998184686; }
+
+size_t Rk4Kennedy::imex_order() const { return 4; }
+
+size_t Rk4Kennedy::implicit_stage_order() const { return 2; }
+
+// The numbers in the tableaus are not exact.  Despite being written
+// as rationals, the true values are irrational numbers.
+const RungeKutta::ButcherTableau& Rk4Kennedy::butcher_tableau() const {
+  static const ButcherTableau tableau{
+      // Substep times
+      {1.0 / 2.0, 83.0 / 250.0, 31.0 / 50.0, 17.0 / 20.0, 1.0},
+      // Substep coefficients
+      {{1.0 / 2.0},
+       {13861.0 / 62500.0, 6889.0 / 62500.0},
+       {-116923316275.0 / 2393684061468.0, -2731218467317.0 / 15368042101831.0,
+        9408046702089.0 / 11113171139209.0},
+       {-451086348788.0 / 2902428689909.0, -2682348792572.0 / 7519795681897.0,
+        12662868775082.0 / 11960479115383.0,
+        3355817975965.0 / 11060851509271.0},
+       {647845179188.0 / 3216320057751.0, 73281519250.0 / 8382639484533.0,
+        552539513391.0 / 3454668386233.0, 3354512671639.0 / 8306763924573.0,
+        4040.0 / 17871.0}},
+      // Result coefficients
+      {82889.0 / 524892.0, 0.0, 15625.0 / 83664.0, 69875.0 / 102672.0,
+       -2260.0 / 8211.0, 1.0 / 4.0},
+      // Coefficients for the embedded method for generating an error measure.
+      {4586570599.0 / 29645900160.0, 0.0, 178811875.0 / 945068544.0,
+       814220225.0 / 1159782912.0, -3700637.0 / 11593932.0, 61727.0 / 225920.0},
+      // Dense output coefficient polynomials
+      {{0.0, 6943876665148.0 / 7220017795957.0, -54480133.0 / 30881146.0,
+        6818779379841.0 / 7100303317025.0},
+       {},
+       {0.0, 7640104374378.0 / 9702883013639.0, -11436875.0 / 14766696.0,
+        2173542590792.0 / 12501825683035.0},
+       {0.0, -20649996744609.0 / 7521556579894.0, 174696575.0 / 18121608.0,
+        -31592104683404.0 / 5083833661969.0},
+       {0.0, 8854892464581.0 / 2390941311638.0, -12120380.0 / 966161.0,
+        61146701046299.0 / 7138195549469.0},
+       {0.0, -11397109935349.0 / 6675773540249.0, 3843.0 / 706.0,
+        -17219254887155.0 / 4939391667607.0}}};
+  return tableau;
+}
+
+const ImexRungeKutta::ImplicitButcherTableau&
+Rk4Kennedy::implicit_butcher_tableau() const {
+  static const ImplicitButcherTableau tableau{
+      {{1.0 / 4.0, 1.0 / 4.0},
+       {8611.0 / 62500.0, -1743.0 / 31250.0, 1.0 / 4.0},
+       {5012029.0 / 34652500.0, -654441.0 / 2922500.0, 174375.0 / 388108.0,
+        1.0 / 4.0},
+       {15267082809.0 / 155376265600.0, -71443401.0 / 120774400.0,
+        730878875.0 / 902184768.0, 2285395.0 / 8070912.0, 1.0 / 4.0},
+       {82889.0 / 524892.0, 0.0, 15625.0 / 83664.0, 69875.0 / 102672.0,
+        -2260.0 / 8211.0, 1.0 / 4.0}}};
+  return tableau;
+}
+}  // namespace TimeSteppers
+
+PUP::able::PUP_ID TimeSteppers::Rk4Kennedy::my_PUP_ID = 0;  // NOLINT

src/Time/TimeSteppers/Rk4Kennedy.hpp

@@ -0,0 +1,64 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "Options/String.hpp"
+#include "Time/TimeSteppers/ImexRungeKutta.hpp"
+#include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace TimeSteppers {
+/*!
+ * \ingroup TimeSteppersGroup
+ * \brief A fourth-order Runge-Kutta method with IMEX support.
+ *
+ * The coefficients are given as ARK4(3)6L[2]SA in \cite Kennedy2003.
+ *
+ * The implicit part is stiffly accurate and L-stable.
+ *
+ * The CFL factor/stable step size is 2.1172491998184686.
+ */
+class Rk4Kennedy : public ImexRungeKutta {
+ public:
+  using options = tmpl::list<>;
+  static constexpr Options::String help = {
+      "A 4th-order Runge-Kutta scheme devised by Kennedy and Carpenter."};
+
+  Rk4Kennedy() = default;
+  Rk4Kennedy(const Rk4Kennedy&) = default;
+  Rk4Kennedy& operator=(const Rk4Kennedy&) = default;
+  Rk4Kennedy(Rk4Kennedy&&) = default;
+  Rk4Kennedy& operator=(Rk4Kennedy&&) = default;
+  ~Rk4Kennedy() override = default;
+
+  size_t order() const override;
+
+  size_t error_estimate_order() const override;
+
+  double stable_step() const override;
+
+  size_t imex_order() const override;
+
+  size_t implicit_stage_order() const override;
+
+  WRAPPED_PUPable_decl_template(Rk4Kennedy);  // NOLINT
+
+  explicit Rk4Kennedy(CkMigrateMessage* /*unused*/) {}
+
+  const ButcherTableau& butcher_tableau() const override;
+
+  const ImplicitButcherTableau& implicit_butcher_tableau() const override;
+};
+
+inline bool constexpr operator==(const Rk4Kennedy& /*lhs*/,
+                                 const Rk4Kennedy& /*rhs*/) {
+  return true;
+}
+
+inline bool constexpr operator!=(const Rk4Kennedy& lhs, const Rk4Kennedy& rhs) {
+  return not(lhs == rhs);
+}
+}  // namespace TimeSteppers

tests/Unit/Helpers/Time/TimeSteppers/RungeKutta.cpp

@@ -18,14 +18,15 @@ namespace {
 // Check order for quadrature (RHS depends only on time).
 void check_quadrature_order(const std::vector<double>& substep_times,
                             const std::vector<double>& coefficients,
+                            const double this_substep_time,
                             const size_t expected) {
   if (expected == 0) {
     return;
   }
   CAPTURE(coefficients);
   CAPTURE(expected);
   // Split out first case to avoid 0^0 annoyance
-  CHECK(alg::accumulate(coefficients, 0.0) == approx(1.0));
+  CHECK(alg::accumulate(coefficients, 0.0) == approx(this_substep_time));
   // Don't require the next order to be inconsistent, as the method
   // may do better for quadrature than for an ODE.  Order 0 (i.e.,
   // that the stepper is at least first order) was checked above.
@@ -36,7 +37,8 @@ void check_quadrature_order(const std::vector<double>& substep_times,
       integral +=
           coefficients[substep] * std::pow(substep_times[substep - 1], order);
     }
-    CHECK(integral == approx(1.0 / (order + 1.0)));
+    CHECK(integral ==
+          approx(pow(this_substep_time, order + 1.0) / (order + 1.0)));
   }
 }
 }  // namespace
@@ -92,8 +94,9 @@ void check_tableau(const TimeSteppers::RungeKutta::ButcherTableau& tableau,
     }
   }
 
-  check_quadrature_order(substep_times, result_coefficients, expected_order);
-  check_quadrature_order(substep_times, error_coefficients,
+  check_quadrature_order(substep_times, result_coefficients, 1.0,
+                         expected_order);
+  check_quadrature_order(substep_times, error_coefficients, 1.0,
                          expected_error_order);
 }
 
@@ -126,7 +129,8 @@ void check_implicit_tableau(
     const auto& coefficients = implicit_coefficients[substep - 1];
     // Substep is DIRK
     CHECK(coefficients.size() <= substep + 1);
-    check_quadrature_order(substep_times, coefficients, expected_stage_order);
+    check_quadrature_order(substep_times, coefficients,
+                           substep_times[substep - 1], expected_stage_order);
   }
 }
 

tests/Unit/Time/TimeSteppers/CMakeLists.txt

@@ -10,7 +10,9 @@ set(LIBRARY_SOURCES
   TimeSteppers/Test_DormandPrince5.cpp
   TimeSteppers/Test_Heun2.cpp
   TimeSteppers/Test_Rk3HesthavenSsp.cpp
+  TimeSteppers/Test_Rk3Kennedy.cpp
   TimeSteppers/Test_Rk3Owren.cpp
+  TimeSteppers/Test_Rk4Kennedy.cpp
   TimeSteppers/Test_Rk4Owren.cpp
   TimeSteppers/Test_Rk5Owren.cpp
   TimeSteppers/Test_Rk5Tsitouras.cpp