Fixup: Problems with interpolator (domain)

And others.
sxs-collaboration · Aug 27, 2024 · a966601 · a966601
1 parent e0c0f81
commit a966601
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 14 deletions.
diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp
@@ -580,7 +580,7 @@ void BinaryWithGravitationalWavesVariables<DataType>::operator()(
     const gsl::not_null<Scalar<DataType>*> rootfinder_bracket_time_upper,
     const gsl::not_null<Cache*> /*cache*/,
     detail::Tags::RootFinderBracketTimeUpper<DataType> /*meta*/) const {
-  get(*rootfinder_bracket_time_upper) = past_time.back();
+  get(*rootfinder_bracket_time_upper) = max_time_interpolator;
 
 }  // namespace detail
 
@@ -901,6 +901,17 @@ void BinaryWithGravitationalWavesVariables<
         static_cast<std::vector<double>>(past_dt_momentum_right.at(i)),
         past_time.front(), std::abs(past_time[0] - past_time[1]));
   }
+  // Get the domain of the interpolation to not trigger domain error on 0
+  // (zero). The maximum time varies by machine roundoff. The interpolation is
+  // done again because above is casted as std::function which does not have
+  // access to the domain.
+  max_time_interpolator =
+      cardinal_cubic_hermite(
+          static_cast<std::vector<double>>(past_position_left.at(2)),
+          static_cast<std::vector<double>>(past_dt_position_left.at(2)),
+          past_time.front(), std::abs(past_time[0] - past_time[1]))
+          .domain()
+          .second;
 }
 
 template <typename DataType>

diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp
@@ -415,6 +415,7 @@ struct BinaryWithGravitationalWavesVariables
                                     const std::array<double, 3>& b) const;
   Scalar<DataType> this_dot_product(const std::array<double, 3>& a,
                                     const tnsr::I<DataType, 3>& b) const;
+  double max_time_interpolator = std::numeric_limits<double>::signaling_NaN();
   void interpolate_past_history();
   DataType find_retarded_time_left(gsl::not_null<Cache*> cache) const;
   DataType find_retarded_time_right(gsl::not_null<Cache*> cache) const;

diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp
@@ -233,19 +233,11 @@ void test_data(const double mass_left, const double mass_right,
     pypp::check_with_random_values<1>(
         &BinaryWithGravitationalWavesProxyExact::test_variables, proxy,
         "BinaryWithGravitationalWaves",
-        {"distance_left",
-         "distance_right",
-         "deriv_one_over_distance_left",
-         "deriv_one_over_distance_right",
-         "deriv_3_distance_left",
-         "deriv_3_distance_right",
-         "normal_left",
-         "normal_right",
-         "near_zone_term",
-         "present_term",
-         "pn_conjugate_momentum3",
-         "pn_extrinsic_curvature",
-         "root_finder_bracket_time_lower",
+        {"distance_left", "distance_right", "deriv_one_over_distance_left",
+         "deriv_one_over_distance_right", "deriv_3_distance_left",
+         "deriv_3_distance_right", "normal_left", "normal_right",
+         "near_zone_term", "present_term", "pn_conjugate_momentum3",
+         "pn_extrinsic_curvature", "root_finder_bracket_time_lower",
          "root_finder_bracket_time_upper",
          "shift_background",
          "longitudinal_shift_background",