Add expiration_after method to FunctionOfTime

Necessary for deterministically coupling the evolution step size to the control system updates.
sxs-collaboration · Oct 11, 2023 · 36a8f0c · 36a8f0c
1 parent b6ed934
commit 36a8f0c
Show file tree

Hide file tree

Showing 12 changed files with 57 additions and 0 deletions.
diff --git a/src/Domain/FunctionsOfTime/FixedSpeedCubic.hpp b/src/Domain/FunctionsOfTime/FixedSpeedCubic.hpp
@@ -70,6 +70,10 @@ class FixedSpeedCubic : public FunctionOfTime {
     return {{initial_time_, std::numeric_limits<double>::infinity()}};
   }
 
+  double expiration_after(const double /*time*/) const override {
+    return std::numeric_limits<double>::infinity();
+  }
+
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override;
 

diff --git a/src/Domain/FunctionsOfTime/FunctionOfTime.hpp b/src/Domain/FunctionsOfTime/FunctionOfTime.hpp
@@ -55,6 +55,16 @@ class FunctionOfTime : public PUP::able {
   /// interval.
   virtual std::array<double, 2> time_bounds() const = 0;
 
+  /// \brief The first expiration time after \p time.
+  ///
+  /// \details For non-updatable functions, this returns infinity. For
+  /// updatable functions, the first expiration time after \p time is
+  /// found by determining the update immediately before \p time. The
+  /// expiration time of this update is what is returned. If \p time
+  /// happens to be an update itself, then the expiration of that
+  /// update is returned.
+  virtual double expiration_after(double time) const = 0;
+
   /// Updates the maximum derivative of the FunctionOfTime at a given time while
   /// also resetting the expiration. By default, a FunctionOfTime cannot be
   /// updated.

diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp
@@ -149,6 +149,12 @@ std::array<double, 2> PiecewisePolynomial<MaxDeriv>::time_bounds() const {
            deriv_info_at_update_times_.expiration_time()}};
 }
 
+template <size_t MaxDeriv>
+double PiecewisePolynomial<MaxDeriv>::expiration_after(
+    const double time) const {
+  return deriv_info_at_update_times_.expiration_after(time);
+}
+
 template <size_t MaxDeriv>
 void PiecewisePolynomial<MaxDeriv>::pup(PUP::er& p) {
   FunctionOfTime::pup(p);

diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp
@@ -78,6 +78,8 @@ class PiecewisePolynomial : public FunctionOfTime {
   /// including the extrapolation region.
   std::array<double, 2> time_bounds() const override;
 
+  double expiration_after(double time) const override;
+
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override;
 

diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp
@@ -83,6 +83,12 @@ std::array<double, 2> QuaternionFunctionOfTime<MaxDeriv>::time_bounds() const {
                     stored_quaternions_and_times_.expiration_time()};
 }
 
+template <size_t MaxDeriv>
+double QuaternionFunctionOfTime<MaxDeriv>::expiration_after(
+    const double time) const {
+  return stored_quaternions_and_times_.expiration_after(time);
+}
+
 template <size_t MaxDeriv>
 void QuaternionFunctionOfTime<MaxDeriv>::pup(PUP::er& p) {
   FunctionOfTime::pup(p);

diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp
@@ -84,6 +84,8 @@ class QuaternionFunctionOfTime : public FunctionOfTime {
   /// Returns domain of validity for the function of time
   std::array<double, 2> time_bounds() const override;
 
+  double expiration_after(double time) const override;
+
   /// Updates the `MaxDeriv`th derivative of the angle piecewisepolynomial at
   /// the given time, then updates the stored quaternions.
   ///

diff --git a/src/Domain/FunctionsOfTime/SettleToConstant.hpp b/src/Domain/FunctionsOfTime/SettleToConstant.hpp
@@ -65,6 +65,10 @@ class SettleToConstant : public FunctionOfTime {
     return {{match_time_, std::numeric_limits<double>::infinity()}};
   }
 
+  double expiration_after(const double /*time*/) const override {
+    return std::numeric_limits<double>::infinity();
+  }
+
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override;
 

diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp
@@ -60,6 +60,8 @@ void test(
   const auto t_bounds = f_of_t->time_bounds();
   CHECK(t_bounds[0] == initial_time);
   CHECK(t_bounds[1] == std::numeric_limits<double>::infinity());
+  CHECK(f_of_t->expiration_after(check_time) ==
+        std::numeric_limits<double>::infinity());
 
   INFO("Test stream operator.");
   CHECK(

diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_OutputTimeBounds.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_OutputTimeBounds.cpp
@@ -42,6 +42,8 @@ class TestFoT : public domain::FunctionsOfTime::FunctionOfTime {
     return {lower_bound_, upper_bound_};
   }
 
+  double expiration_after(const double /*time*/) const override { ERROR(""); }
+
   std::array<DataVector, 1> func(const double /*t*/) const override {
     ERROR("");
   }

diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp
@@ -51,12 +51,16 @@ void test(const gsl::not_null<FunctionsOfTime::FunctionOfTime*> f_of_t,
 
     t += dt;
     f_of_t_derived->update(t, {6.0, 0.0}, t + dt);
+    CHECK(f_of_t->expiration_after(t) == t + dt);
+    CHECK(f_of_t->expiration_after(t + 0.5 * dt) == t + dt);
+    CHECK(f_of_t->expiration_after(t - 0.5 * dt) == t);
     CHECK(*f_of_t_derived != f_of_t_derived_copy);
   }
   // test time_bounds function
   const auto t_bounds = f_of_t->time_bounds();
   CHECK(t_bounds[0] == 0.0);
   CHECK(t_bounds[1] == 4.8);
+  CHECK(f_of_t->expiration_after(0.0) == dt);
 }
 
 template <size_t DerivOrder>

diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp
@@ -280,7 +280,11 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime",
       t += time_step;
       expir_time += time_step;
       qfot.update(t, DataVector{3, 0.0}, expir_time);
+      CHECK(qfot.expiration_after(t) == expir_time);
+      CHECK(qfot.expiration_after(t + 0.5 * time_step) == expir_time);
+      CHECK(qfot.expiration_after(t - 0.5 * time_step) == t);
     }
+    CHECK(qfot.expiration_after(0.0) == time_step);
 
     // Get the quaternion and 2 derivatives at a certain time.
     double check_time = 5.398;
@@ -352,7 +356,11 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime",
       t += time_step;
       expir_time += time_step;
       qfot.update(t, DataVector{{0.0, 0.0, 2.0 * fac1}}, expir_time);
+      CHECK(qfot.expiration_after(t) == expir_time);
+      CHECK(qfot.expiration_after(t + 0.5 * time_step) == expir_time);
+      CHECK(qfot.expiration_after(t - 0.5 * time_step) == t);
     }
+    CHECK(qfot.expiration_after(0.0) == time_step);
 
     // Get the quaternion and 2 derivatives at a certain time.
     double check_time = 5.398;
@@ -431,7 +439,11 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime",
       t += time_step;
       expir_time += time_step;
       qfot.update(t, DataVector{{0.0, 0.0, 6.0 * fac1}}, expir_time);
+      CHECK(qfot.expiration_after(t) == expir_time);
+      CHECK(qfot.expiration_after(t + 0.5 * time_step) == expir_time);
+      CHECK(qfot.expiration_after(t - 0.5 * time_step) == t);
     }
+    CHECK(qfot.expiration_after(0.0) == time_step);
 
     // Get the quaternion and 2 derivatives at a certain time.
     double check_time = 5.398;

diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp
@@ -47,6 +47,9 @@ void test(
   // test time_bounds function
   const auto t_bounds = f_of_t->time_bounds();
   CHECK(t_bounds[0] == 10.0);
+  CHECK(t_bounds[1] == std::numeric_limits<double>::infinity());
+  CHECK(f_of_t->expiration_after(match_time) ==
+        std::numeric_limits<double>::infinity());
 }
 }  // namespace