Change Heun2 dense output to linear interpolation

The quadratic output was highly oscillatory in IMEX evolutions.
sxs-collaboration · Feb 5, 2024 · 77dc0bc · 77dc0bc
1 parent a5f4edd
commit 77dc0bc
Show file tree

Hide file tree

Showing 17 changed files with 75 additions and 24 deletions.
diff --git a/src/Time/TimeSteppers/Heun2.cpp b/src/Time/TimeSteppers/Heun2.cpp
@@ -31,8 +31,8 @@ const RungeKutta::ButcherTableau& Heun2::butcher_tableau() const {
       // Coefficients for the embedded method for generating an error measure.
       {1.0, 0.0},
       // Dense output coefficient polynomials
-      {{0.0, 1.0, -0.5},
-       {0.0, 0.0, 0.5}}};
+      {{0.0, 0.5},
+       {0.0, 0.5}}};
   return tableau;
 }
 

diff --git a/tests/Unit/Helpers/Time/TimeSteppers/ImexHelpers.cpp b/tests/Unit/Helpers/Time/TimeSteppers/ImexHelpers.cpp
@@ -86,4 +86,49 @@ void check_convergence_order(const ImexTimeStepper& stepper,
   CHECK(convergence_rate(step_range, 1, do_integral, output) ==
         approx(stepper.imex_order()).margin(0.4));
 }
+
+void check_bounded_dense_output(const ImexTimeStepper& stepper) {
+  const double decay_constant = 1.0e5;
+  const Slab slab(0.0, 1.0);
+  const auto time_step = slab.duration();
+  const auto rhs = [&](const double v, const double /*t*/) {
+    return -decay_constant * v;
+  };
+
+  TimeStepId time_step_id(true, 0, slab.start());
+  double y = 1.0;
+  TimeSteppers::History<double> history{stepper.order()};
+  initialize_history(
+      time_step_id.step_time(), make_not_null(&history),
+      [&](const double t) { return exp(-decay_constant * t); }, rhs, time_step,
+      stepper.number_of_past_steps());
+  double test_time = 0.0;
+  for (;;) {
+    history.insert(time_step_id, TimeSteppers::History<double>::no_value,
+                   rhs(y, time_step_id.substep_time()));
+
+    while (test_time < 1.0) {
+      y = 1.0;
+      if (not stepper.dense_update_u(make_not_null(&y), history, test_time)) {
+        break;
+      }
+      CHECK(abs(y) < 5.0);
+      test_time += 0.1;
+    }
+
+    if (time_step_id.slab_number() != 0) {
+      break;
+    }
+
+    y = 1.0;
+    stepper.add_inhomogeneous_implicit_terms(
+        make_not_null(&y), make_not_null(&history), time_step);
+    const double weight =
+        stepper.implicit_weight(make_not_null(&history), time_step);
+
+    y /= 1.0 + decay_constant * weight;
+    time_step_id = stepper.next_time_id(time_step_id, time_step);
+  }
+  CHECK(test_time >= 1.0);
+}
 }  // namespace TimeStepperTestUtils::imex
diff --git a/tests/Unit/Helpers/Time/TimeSteppers/ImexHelpers.hpp b/tests/Unit/Helpers/Time/TimeSteppers/ImexHelpers.hpp
@@ -14,4 +14,7 @@ namespace TimeStepperTestUtils::imex {
 void check_convergence_order(const ImexTimeStepper& stepper,
                              const std::pair<int32_t, int32_t>& step_range,
                              bool output = false);
+
+/// Check that dense output does not have large oscillations.
+void check_bounded_dense_output(const ImexTimeStepper& stepper);
 }  // namespace TimeStepperTestUtils::imex
diff --git a/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp b/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp
@@ -339,8 +339,9 @@ void check_convergence_order(const TimeStepper& stepper,
         approx(stepper.order()).margin(0.4));
 }
 
-void check_dense_output(const TimeStepper& stepper,
-                        const size_t history_integration_order) {
+void check_dense_output(
+    const TimeStepper& stepper, const size_t history_integration_order,
+    const std::pair<int32_t, int32_t>& convergence_step_range) {
   const auto get_dense = [&stepper, &history_integration_order](
                              const TimeDelta& step_size, const double time) {
     const auto impl = [&stepper, &history_integration_order, &step_size,
@@ -418,24 +419,20 @@ void check_dense_output(const TimeStepper& stepper,
 
   // Test convergence
   {
-    const int32_t large_steps = 10;
-    // The high-order solvers have round-off error around here
-    const int32_t small_steps = 30;
-
     const auto error = [&get_dense](const int32_t steps) {
       const Slab slab(0., 1.);
       return abs(get_dense(slab.duration() / steps, 0.25 * M_PI) -
                  exp(0.25 * M_PI));
     };
-    CHECK(convergence_rate({large_steps, small_steps}, 1, error) ==
+    CHECK(convergence_rate(convergence_step_range, 1, error) ==
           approx(history_integration_order).margin(0.4));
 
     const auto error_backwards = [&get_dense](const int32_t steps) {
       const Slab slab(-1., 0.);
       return abs(get_dense(-slab.duration() / steps, -0.25 * M_PI) -
                  exp(-0.25 * M_PI));
     };
-    CHECK(convergence_rate({large_steps, small_steps}, 1, error_backwards) ==
+    CHECK(convergence_rate(convergence_step_range, 1, error_backwards) ==
           approx(history_integration_order).margin(0.4));
   }
 }

diff --git a/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.hpp b/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.hpp
@@ -91,8 +91,9 @@ void check_convergence_order(const TimeStepper& stepper,
                              const std::pair<int32_t, int32_t>& step_range,
                              bool output = false);
 
-void check_dense_output(const TimeStepper& stepper,
-                        const size_t history_integration_order);
+void check_dense_output(
+    const TimeStepper& stepper, const size_t history_integration_order,
+    const std::pair<int32_t, int32_t>& convergence_step_range);
 
 void check_strong_stability_preservation(const TimeStepper& stepper,
                                          double step_size);

diff --git a/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp b/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp
@@ -55,7 +55,8 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth", "[Unit][Time]") {
     TimeStepperTestUtils::check_convergence_order(stepper, {10, 30});
     for(size_t history_order = 1; history_order <= order; ++history_order){
       CAPTURE(history_order);
-      TimeStepperTestUtils::check_dense_output(stepper, history_order);
+      TimeStepperTestUtils::check_dense_output(stepper, history_order,
+                                               {10, 30});
     }
 
     CHECK(stepper.order() == order);

diff --git a/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp b/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp
@@ -56,7 +56,8 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsMoultonPc", "[Unit][Time]") {
     TimeStepperTestUtils::check_convergence_order(stepper, {10, 30});
     for (size_t history_order = 2; history_order <= order; ++history_order) {
       CAPTURE(history_order);
-      TimeStepperTestUtils::check_dense_output(stepper, history_order);
+      TimeStepperTestUtils::check_dense_output(stepper, history_order,
+                                               {10, 30});
     }
   }
 

diff --git a/tests/Unit/Time/TimeSteppers/Test_ClassicalRungeKutta4.cpp b/tests/Unit/Time/TimeSteppers/Test_ClassicalRungeKutta4.cpp
@@ -33,7 +33,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.ClassicalRungeKutta4",
   TimeStepperTestUtils::integrate_variable_test(stepper, 4, 0, 1.0e-9);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
   TimeStepperTestUtils::stability_test(stepper);
-  TimeStepperTestUtils::check_dense_output(stepper, 4_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 4_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper,
                                      TimeSteppers::ClassicalRungeKutta4>(

diff --git a/tests/Unit/Time/TimeSteppers/Test_DormandPrince5.cpp b/tests/Unit/Time/TimeSteppers/Test_DormandPrince5.cpp
@@ -32,7 +32,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.DormandPrince5", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 5, 0, 1.0e-9);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
   TimeStepperTestUtils::stability_test(stepper);
-  TimeStepperTestUtils::check_dense_output(stepper, 5_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 5_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::DormandPrince5>(
       "DormandPrince5");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Heun2.cpp b/tests/Unit/Time/TimeSteppers/Test_Heun2.cpp
@@ -33,9 +33,10 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Heun2", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 2, 0, 1.0e-6);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 2_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 2_st, {10, 100});
 
   TimeStepperTestUtils::imex::check_convergence_order(stepper, {10, 50});
+  TimeStepperTestUtils::imex::check_bounded_dense_output(stepper);
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Heun2>("Heun2");
   test_serialization(stepper);

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk3HesthavenSsp.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk3HesthavenSsp.cpp
@@ -31,7 +31,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk3HesthavenSsp", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 3, 0, 1e-9);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 3_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 3_st, {10, 30});
 
   TimeStepperTestUtils::check_strong_stability_preservation(stepper, 1.0);
 

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk3Kennedy.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk3Kennedy.cpp
@@ -33,9 +33,10 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk3Kennedy", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 3, 0, 1.0e-9);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 3_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 3_st, {10, 30});
 
   TimeStepperTestUtils::imex::check_convergence_order(stepper, {10, 50});
+  TimeStepperTestUtils::imex::check_bounded_dense_output(stepper);
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk3Kennedy>(
       "Rk3Kennedy");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk3Owren.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk3Owren.cpp
@@ -31,7 +31,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk3Owren", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 3, 0, 1.0e-9);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 3_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 3_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk3Owren>(
       "Rk3Owren");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk4Kennedy.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk4Kennedy.cpp
@@ -33,9 +33,10 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk4Kennedy", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 4, 0, 1.0e-14);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 4_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 4_st, {10, 30});
 
   TimeStepperTestUtils::imex::check_convergence_order(stepper, {10, 50});
+  TimeStepperTestUtils::imex::check_bounded_dense_output(stepper);
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk4Kennedy>(
       "Rk4Kennedy");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk4Owren.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk4Owren.cpp
@@ -31,7 +31,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk4Owren", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 4, 0, 1.0e-9);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 4_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 4_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk4Owren>(
       "Rk4Owren");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk5Owren.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk5Owren.cpp
@@ -31,7 +31,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk5Owren", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 5, 0, 1.0e-9);
   TimeStepperTestUtils::stability_test(stepper);
   TimeStepperTestUtils::check_convergence_order(stepper, {10, 50});
-  TimeStepperTestUtils::check_dense_output(stepper, 5_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 5_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk5Owren>(
       "Rk5Owren");

diff --git a/tests/Unit/Time/TimeSteppers/Test_Rk5Tsitouras.cpp b/tests/Unit/Time/TimeSteppers/Test_Rk5Tsitouras.cpp
@@ -32,7 +32,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Rk5Tsitouras", "[Unit][Time]") {
   TimeStepperTestUtils::integrate_variable_test(stepper, 5, 0, 1.0e-9);
   TimeStepperTestUtils::check_convergence_order(stepper, {30, 70});
   TimeStepperTestUtils::stability_test(stepper);
-  TimeStepperTestUtils::check_dense_output(stepper, 5_st);
+  TimeStepperTestUtils::check_dense_output(stepper, 5_st, {10, 30});
 
   TestHelpers::test_factory_creation<TimeStepper, TimeSteppers::Rk5Tsitouras>(
       "Rk5Tsitouras");