SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/Cce - ScriPlusInterpolationManager.hpp Hit Total Coverage
Commit: 817e13c5144619b701c7cd870655d8dbf94ab8ce Lines: 24 45 53.3 %
Date: 2024-07-19 22:17:05
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <complex>
       8             : #include <cstddef>
       9             : #include <deque>
      10             : #include <memory>
      11             : #include <pup.h>
      12             : #include <pup_stl.h>
      13             : #include <utility>
      14             : #include <vector>
      15             : 
      16             : #include "DataStructures/ApplyMatrices.hpp"
      17             : #include "DataStructures/DataVector.hpp"
      18             : #include "Evolution/Systems/Cce/WorldtubeDataManager.hpp"
      19             : #include "NumericalAlgorithms/Interpolation/SpanInterpolator.hpp"
      20             : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      21             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      22             : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
      23             : #include "Utilities/Algorithm.hpp"
      24             : #include "Utilities/ErrorHandling/Assert.hpp"
      25             : #include "Utilities/ErrorHandling/Error.hpp"
      26             : #include "Utilities/Literals.hpp"
      27             : #include "Utilities/MakeArray.hpp"
      28             : 
      29             : namespace Cce {
      30             : 
      31             : /*!
      32             :  * \brief Stores necessary data and interpolates on to new time points at scri+.
      33             :  *
      34             :  * \details The coordinate time used for the CCE evolution is not the same as
      35             :  * the asymptotic inertial retarded time, which is determined through a separate
      36             :  * evolution equation. This class manages the scri+ data passed in (via
      37             :  * `insert_data()`) along with the set of inertial retarded times associated
      38             :  * with that data, and interpolates to a set of requested times (supplied via
      39             :  * `insert_target_time()`).
      40             :  *
      41             :  * Template parameters:
      42             :  * - `VectorTypeToInterpolate`: the vector type associated with the values to
      43             :  * interpolate.
      44             :  * - `Tag`: The tag associated with the interpolation procedure. This determines
      45             :  * the behavior of the interpolation return value. The default is just
      46             :  * interpolation, if `Tag` has prefix `::Tags::Multiplies` or `Tags::Du`, the
      47             :  * interpolator performs the additional multiplication or time derivative as a
      48             :  * step in the interpolation procedure.
      49             :  */
      50             : template <typename VectorTypeToInterpolate, typename Tag>
      51           1 : struct ScriPlusInterpolationManager {
      52             :  public:
      53           0 :   ScriPlusInterpolationManager() = default;
      54             : 
      55           0 :   ScriPlusInterpolationManager(
      56             :       const size_t target_number_of_points, const size_t vector_size,
      57             :       std::unique_ptr<intrp::SpanInterpolator> interpolator)
      58             :       : vector_size_{vector_size},
      59             :         target_number_of_points_{target_number_of_points},
      60             :         interpolator_{std::move(interpolator)} {}
      61             : 
      62             :   /// \brief provide data to the interpolation manager.
      63             :   ///
      64             :   /// \details `u_bondi` is a vector of inertial times, and `to_interpolate` is
      65             :   /// the vector of values that will be interpolated to target times.
      66           1 :   void insert_data(const DataVector& u_bondi,
      67             :                    const VectorTypeToInterpolate& to_interpolate) {
      68             :     ASSERT(to_interpolate.size() == vector_size_,
      69             :            "Inserted data must be of size specified at construction: "
      70             :                << vector_size_
      71             :                << " and provided data is of size: " << to_interpolate.size());
      72             :     u_bondi_values_.push_back(u_bondi);
      73             :     to_interpolate_values_.push_back(to_interpolate);
      74             :     u_bondi_ranges_.emplace_back(min(u_bondi), max(u_bondi));
      75             :   }
      76             : 
      77             :   /// \brief Request a target time to be interpolated to when enough data has
      78             :   /// been accumulated.
      79             :   ///
      80             :   /// For optimization, we assume that these are inserted in ascending order.
      81           1 :   void insert_target_time(const double time) {
      82             :     target_times_.push_back(time);
      83             :     // if this is the first time in the deque, we should use it to determine
      84             :     // whether some of the old data is ready to be removed. In other cases, this
      85             :     // check can be performed when a time is removed from the queue via
      86             :     // `interpolate_and_pop_first_time`
      87             :     if (target_times_.size() == 1) {
      88             :       remove_unneeded_early_times();
      89             :     }
      90             :   }
      91             : 
      92             :   /// \brief Determines whether enough data before and after the first time in
      93             :   /// the target time queue has been provided to interpolate.
      94             :   ///
      95             :   /// \details If possible, this function will require that the target time to
      96             :   /// be interpolated is reasonably centered on the range, but will settle for
      97             :   /// non-centered data if the time is too early for the given data, which is
      98             :   /// necessary to get the first couple of times out of the simulation. This
      99             :   /// function always returns false if all of the provided data is earlier than
     100             :   /// the next target time, indicating that the caller should wait until more
     101             :   /// data has been provided before interpolating.
     102           1 :   bool first_time_is_ready_to_interpolate() const;
     103             : 
     104           0 :   const std::deque<std::pair<double, double>>& get_u_bondi_ranges() const {
     105             :     return u_bondi_ranges_;
     106             :   }
     107             : 
     108             :   /// \brief Interpolate to the first target time in the queue, returning both
     109             :   /// the time and the interpolated data at that time.
     110             :   ///
     111             :   /// \note If this function is not able to interpolate to the full accuracy
     112             :   /// using a centered stencil, it will perform the best interpolation
     113             :   /// available.
     114             :   ///
     115             :   /// \warning This function will extrapolate if the target times are
     116             :   /// outside the range of data the class has been provided. This is intentional
     117             :   /// to support small extrapolation at the end of a simulation when no further
     118             :   /// data is available, but for full accuracy, check
     119             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     120             :   /// functions
     121           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_first_time();
     122             : 
     123             :   /// \brief Interpolate to the first target time in the queue, returning both
     124             :   /// the time and the interpolated data at that time, and remove the first time
     125             :   /// from the queue.
     126             :   ///
     127             :   /// \note If this function is not able to interpolate to the full accuracy
     128             :   /// using a centered stencil, it will perform the best interpolation
     129             :   /// available.
     130             :   ///
     131             :   /// \warning This function will extrapolate if the target times are
     132             :   /// outside the range of data the class has been provided. This is intentional
     133             :   /// to support small extrapolation at the end of a simulation when no further
     134             :   /// data is available, but for full accuracy, check
     135             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     136             :   /// functions
     137           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_and_pop_first_time();
     138             : 
     139             :   /// \brief return the number of times in the target times queue
     140           1 :   size_t number_of_target_times() const { return target_times_.size(); }
     141             : 
     142             :   /// \brief return the number of data points that have been provided to the
     143             :   /// interpolation manager
     144           1 :   size_t number_of_data_points() const { return u_bondi_ranges_.size(); }
     145             : 
     146             :   // NOLINTNEXTLINE(google-runtime-references)
     147           0 :   void pup(PUP::er& p) {
     148             :     p | u_bondi_values_;
     149             :     p | to_interpolate_values_;
     150             :     p | u_bondi_ranges_;
     151             :     p | target_times_;
     152             :     p | vector_size_;
     153             :     p | target_number_of_points_;
     154             :     p | interpolator_;
     155             :   }
     156             : 
     157             :  private:
     158           0 :   void remove_unneeded_early_times();
     159             : 
     160             :   friend struct ScriPlusInterpolationManager<VectorTypeToInterpolate,
     161             :                                              Tags::Du<Tag>>;
     162             : 
     163           0 :   std::deque<DataVector> u_bondi_values_;
     164           0 :   std::deque<VectorTypeToInterpolate> to_interpolate_values_;
     165           0 :   std::deque<std::pair<double, double>> u_bondi_ranges_;
     166           0 :   std::deque<double> target_times_;
     167           0 :   size_t vector_size_ = 0_st;
     168           0 :   size_t target_number_of_points_ = 0_st;
     169           0 :   std::unique_ptr<intrp::SpanInterpolator> interpolator_;
     170             : };
     171             : 
     172             : template <typename VectorTypeToInterpolate, typename Tag>
     173             : bool ScriPlusInterpolationManager<
     174             :     VectorTypeToInterpolate, Tag>::first_time_is_ready_to_interpolate() const {
     175             :   if (target_times_.empty()) {
     176             :     return false;
     177             :   }
     178             :   auto maxes_below = alg::count_if(
     179             :       u_bondi_ranges_, [this](const std::pair<double, double> time) {
     180             :         return time.second <= target_times_.front();
     181             :       });
     182             :   auto mins_above = alg::count_if(u_bondi_ranges_,
     183             :                                   [this](const std::pair<double, double> time) {
     184             :                                     return time.first > target_times_.front();
     185             :                                   });
     186             : 
     187             :   // we might ask for a time that's too close to the end or the beginning of
     188             :   // our data, in which case we will settle for at least one point below and
     189             :   // above and a sufficient number of total points.
     190             :   // This will always be `false` if the earliest target time is later than all
     191             :   // provided data points, which would require extrapolation
     192             :   return ((static_cast<size_t>(maxes_below) > target_number_of_points_ or
     193             :            static_cast<size_t>(maxes_below + mins_above) >
     194             :                2 * target_number_of_points_) and
     195             :           static_cast<size_t>(mins_above) > target_number_of_points_);
     196             : }
     197             : 
     198             : template <typename VectorTypeToInterpolate, typename Tag>
     199             : std::pair<double, VectorTypeToInterpolate> ScriPlusInterpolationManager<
     200             :     VectorTypeToInterpolate, Tag>::interpolate_first_time() {
     201             :   if (target_times_.empty()) {
     202             :     ERROR("There are no target times to interpolate.");
     203             :   }
     204             :   if (to_interpolate_values_.size() < 2 * target_number_of_points_) {
     205             :     ERROR("Insufficient data points to continue interpolation: have "
     206             :           << to_interpolate_values_.size() << ", need at least"
     207             :           << 2 * target_number_of_points_);
     208             :   }
     209             : 
     210             :   VectorTypeToInterpolate result{vector_size_};
     211             :   const size_t interpolation_data_size = to_interpolate_values_.size();
     212             : 
     213             :   VectorTypeToInterpolate interpolation_values{2 * target_number_of_points_};
     214             :   DataVector interpolation_times{2 * target_number_of_points_};
     215             :   for (size_t i = 0; i < vector_size_; ++i) {
     216             :     // binary search assumes times placed in sorted order
     217             :     auto upper_bound_offset = static_cast<size_t>(std::distance(
     218             :         u_bondi_values_.begin(),
     219             :         std::upper_bound(u_bondi_values_.begin(), u_bondi_values_.end(),
     220             :                          target_times_.front(),
     221             :                          [&i](const double rhs, const DataVector& lhs) {
     222             :                            return rhs < lhs[i];
     223             :                          })));
     224             :     size_t lower_bound_offset =
     225             :         upper_bound_offset == 0 ? 0 : upper_bound_offset - 1;
     226             : 
     227             :     if (upper_bound_offset + target_number_of_points_ >
     228             :         interpolation_data_size) {
     229             :       upper_bound_offset = interpolation_data_size;
     230             :       lower_bound_offset =
     231             :           interpolation_data_size - 2 * target_number_of_points_;
     232             :     } else if (lower_bound_offset < target_number_of_points_ - 1) {
     233             :       lower_bound_offset = 0;
     234             :       upper_bound_offset = 2 * target_number_of_points_;
     235             :     } else {
     236             :       lower_bound_offset = lower_bound_offset + 1 - target_number_of_points_;
     237             :       upper_bound_offset = lower_bound_offset + 2 * target_number_of_points_;
     238             :     }
     239             :     auto interpolation_values_begin =
     240             :         to_interpolate_values_.begin() +
     241             :         static_cast<ptrdiff_t>(lower_bound_offset);
     242             :     auto interpolation_times_begin =
     243             :         u_bondi_values_.begin() + static_cast<ptrdiff_t>(lower_bound_offset);
     244             :     auto interpolation_times_end =
     245             :         u_bondi_values_.begin() + static_cast<ptrdiff_t>(upper_bound_offset);
     246             : 
     247             :     // interpolate using the data sets in the restricted iterators
     248             :     auto value_it = interpolation_values_begin;
     249             :     size_t vector_position = 0;
     250             :     for (auto time_it = interpolation_times_begin;
     251             :          time_it != interpolation_times_end;
     252             :          ++time_it, ++value_it, ++vector_position) {
     253             :       interpolation_values[vector_position] = (*value_it)[i];
     254             :       interpolation_times[vector_position] = (*time_it)[i];
     255             :     }
     256             :     result[i] = interpolator_->interpolate(
     257             :         gsl::span<const double>(interpolation_times.data(),
     258             :                                 interpolation_times.size()),
     259             :         gsl::span<const typename VectorTypeToInterpolate::value_type>(
     260             :             interpolation_values.data(), interpolation_values.size()),
     261             :         target_times_.front());
     262             :   }
     263             :   return std::make_pair(target_times_.front(), std::move(result));
     264             : }
     265             : 
     266             : template <typename VectorTypeToInterpolate, typename Tag>
     267             : std::pair<double, VectorTypeToInterpolate> ScriPlusInterpolationManager<
     268             :     VectorTypeToInterpolate, Tag>::interpolate_and_pop_first_time() {
     269             :   std::pair<double, VectorTypeToInterpolate> interpolated =
     270             :       interpolate_first_time();
     271             :   target_times_.pop_front();
     272             : 
     273             :   if (not target_times_.empty()) {
     274             :     remove_unneeded_early_times();
     275             :   }
     276             :   return interpolated;
     277             : }
     278             : 
     279             : template <typename VectorTypeToInterpolate, typename Tag>
     280             : void ScriPlusInterpolationManager<VectorTypeToInterpolate,
     281             :                                   Tag>::remove_unneeded_early_times() {
     282             :   // pop times we no longer need because their maxes are too far in the past
     283             :   auto time_it = u_bondi_ranges_.begin();
     284             :   size_t times_counter = 0;
     285             :   while (time_it < u_bondi_ranges_.end() and
     286             :          (*time_it).second < target_times_.front()) {
     287             :     if (times_counter > target_number_of_points_ and
     288             :         u_bondi_ranges_.size() >= 2 * target_number_of_points_) {
     289             :       u_bondi_ranges_.pop_front();
     290             :       u_bondi_values_.pop_front();
     291             :       to_interpolate_values_.pop_front();
     292             :     } else {
     293             :       ++times_counter;
     294             :     }
     295             :     ++time_it;
     296             :   }
     297             : }
     298             : 
     299             : /*!
     300             :  * \brief Stores necessary data and interpolates on to new time points at scri+,
     301             :  * multiplying two results together before supplying the result.
     302             :  *
     303             :  * \details The coordinate time used for the CCE evolution is not the same as
     304             :  * the asymptotic inertial retarded time, which is determined through a separate
     305             :  * evolution equation. This class manages the two sets of  scri+ data passed in
     306             :  * (via `insert_data()`) along with the set of inertial retarded times
     307             :  * associated with that data, and interpolates to a set of inertial requested
     308             :  * times (supplied via `insert_target_time()`), multiplying the two
     309             :  * interpolation results together before returning.
     310             :  *
     311             :  * Template parameters:
     312             :  * - `VectorTypeToInterpolate`: the vector type associated with the values to
     313             :  * interpolate.
     314             :  * - `::Tags::Multiplies<MultipliesLhs, MultipliesRhs>`: The tag associated with
     315             :  * the interpolation procedure. This determines the behavior of the
     316             :  * interpolation return value. The default is just interpolation, if `Tag` has
     317             :  * prefix  `::Tags::Multiplies` (this case) or `Tags::Du`, the interpolator
     318             :  * performs the additional multiplication or time derivative as a step in the
     319             :  * interpolation procedure.
     320             :  */
     321             : template <typename VectorTypeToInterpolate, typename MultipliesLhs,
     322             :           typename MultipliesRhs>
     323           1 : struct ScriPlusInterpolationManager<
     324             :     VectorTypeToInterpolate, ::Tags::Multiplies<MultipliesLhs, MultipliesRhs>> {
     325             :  public:
     326           0 :   ScriPlusInterpolationManager() = default;
     327           0 :   ScriPlusInterpolationManager(
     328             :       const size_t target_number_of_points, const size_t vector_size,
     329             :       std::unique_ptr<intrp::SpanInterpolator> interpolator)
     330             :       : interpolation_manager_lhs_{target_number_of_points, vector_size,
     331             :                                    interpolator->get_clone()},
     332             :         interpolation_manager_rhs_{target_number_of_points, vector_size,
     333             :                                    std::move(interpolator)} {}
     334             : 
     335             :   /// \brief provide data to the interpolation manager.
     336             :   ///
     337             :   /// \details `u_bondi` is a vector of inertial times, and `to_interpolate_lhs`
     338             :   /// and `to_interpolate_rhs` are the vector of values that will be
     339             :   /// interpolated to target times. The interpolation result will be the product
     340             :   /// of the interpolated lhs and rhs vectors.
     341           1 :   void insert_data(const DataVector& u_bondi,
     342             :                    const VectorTypeToInterpolate& to_interpolate_lhs,
     343             :                    const VectorTypeToInterpolate& to_interpolate_rhs) {
     344             :     interpolation_manager_lhs_.insert_data(u_bondi, to_interpolate_lhs);
     345             :     interpolation_manager_rhs_.insert_data(u_bondi, to_interpolate_rhs);
     346             :   }
     347             : 
     348             :   /// \brief Request a target time to be interpolated to when enough data has
     349             :   /// been accumulated.
     350             :   ///
     351             :   /// For optimization, we assume that these are inserted in ascending order.
     352           1 :   void insert_target_time(const double time) {
     353             :     interpolation_manager_lhs_.insert_target_time(time);
     354             :     interpolation_manager_rhs_.insert_target_time(time);
     355             :   }
     356             : 
     357             :   /// \brief Determines whether enough data before and after the first time in
     358             :   /// the target time queue has been provided to interpolate.
     359             :   ///
     360             :   /// \details If possible, this function will require that the target time to
     361             :   /// be interpolated is reasonably centered on the range, but will settle for
     362             :   /// non-centered data if the time is too early for the given data, which is
     363             :   /// necessary to get the first couple of times out of the simulation. This
     364             :   /// function always returns false if all of the provided data is earlier than
     365             :   /// the next target time, indicating that the caller should wait until more
     366             :   /// data has been provided before interpolating.
     367           1 :   bool first_time_is_ready_to_interpolate() const {
     368             :     return interpolation_manager_lhs_.first_time_is_ready_to_interpolate() and
     369             :            interpolation_manager_rhs_.first_time_is_ready_to_interpolate();
     370             :   }
     371             : 
     372           0 :   const std::deque<std::pair<double, double>>& get_u_bondi_ranges() const {
     373             :     return interpolation_manager_lhs_.get_u_bondi_ranges();
     374             :   }
     375             : 
     376             :   /// \brief Interpolate to the first target time in the queue, returning both
     377             :   /// the time and the interpolated-and-multiplied data at that time.
     378             :   ///
     379             :   /// \note If this function is not able to interpolate to the full accuracy
     380             :   /// using a centered stencil, it will perform the best interpolation
     381             :   /// available.
     382             :   ///
     383             :   /// \warning This function will extrapolate if the target times are
     384             :   /// outside the range of data the class has been provided. This is intentional
     385             :   /// to support small extrapolation at the end of a simulation when no further
     386             :   /// data is available, but for full accuracy, check
     387             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     388             :   /// functions
     389           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_first_time() {
     390             :     const auto lhs_interpolation =
     391             :         interpolation_manager_lhs_.interpolate_first_time();
     392             :     const auto rhs_interpolation =
     393             :         interpolation_manager_rhs_.interpolate_first_time();
     394             :     return std::make_pair(lhs_interpolation.first,
     395             :                           lhs_interpolation.second * rhs_interpolation.second);
     396             :   }
     397             : 
     398             :   /// \brief Interpolate to the first target time in the queue, returning both
     399             :   /// the time and the interpolated-and-multiplied data at that time, and remove
     400             :   /// the first time from the queue.
     401             :   ///
     402             :   /// \note If this function is not able to interpolate to the full accuracy
     403             :   /// using a centered stencil, it will perform the best interpolation
     404             :   /// available.
     405             :   ///
     406             :   /// \warning This function will extrapolate if the target times are
     407             :   /// outside the range of data the class has been provided. This is intentional
     408             :   /// to support small extrapolation at the end of a simulation when no further
     409             :   /// data is available, but for full accuracy, check
     410             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     411             :   /// functions
     412           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_and_pop_first_time() {
     413             :     const auto lhs_interpolation =
     414             :         interpolation_manager_lhs_.interpolate_and_pop_first_time();
     415             :     const auto rhs_interpolation =
     416             :         interpolation_manager_rhs_.interpolate_and_pop_first_time();
     417             :     return std::make_pair(lhs_interpolation.first,
     418             :                           lhs_interpolation.second * rhs_interpolation.second);
     419             :   }
     420             : 
     421             :   /// \brief return the number of times in the target times queue
     422           1 :   size_t number_of_target_times() const {
     423             :     return interpolation_manager_lhs_.number_of_target_times();
     424             :   }
     425             : 
     426             :   /// \brief return the number of data points that have been provided to the
     427             :   /// interpolation manager
     428           1 :   size_t number_of_data_points() const {
     429             :     return interpolation_manager_lhs_.number_of_data_points();
     430             :   }
     431             : 
     432             :   // NOLINTNEXTLINE(google-runtime-references)
     433           0 :   void pup(PUP::er& p) {
     434             :     p | interpolation_manager_lhs_;
     435             :     p | interpolation_manager_rhs_;
     436             :   }
     437             : 
     438             :  private:
     439             :   /// \cond
     440             :   ScriPlusInterpolationManager<VectorTypeToInterpolate, MultipliesLhs>
     441             :       interpolation_manager_lhs_;
     442             :   ScriPlusInterpolationManager<VectorTypeToInterpolate, MultipliesRhs>
     443             :       interpolation_manager_rhs_;
     444             :   /// \endcond
     445             : };
     446             : 
     447             : /*!
     448             :  * \brief Stores necessary data and interpolates on to new time points at scri+,
     449             :  * differentiating before supplying the result.
     450             :  *
     451             :  * \details The coordinate time used for the CCE evolution is not the same as
     452             :  * the asymptotic inertial retarded time, which is determined through a separate
     453             :  * evolution equation. This class manages the scri+ data passed in
     454             :  * (via `insert_data()`) along with the set of inertial retarded times
     455             :  * associated with that data, and interpolates to a set of requested times
     456             :  * (supplied via `insert_target_time()`), differentiating the interpolation
     457             :  * before returning.
     458             :  *
     459             :  * Template parameters:
     460             :  * - `VectorTypeToInterpolate`: the vector type associated with the values to
     461             :  * interpolate.
     462             :  * - `Tags::Du<Tag>`: The tag associated with the interpolation procedure. This
     463             :  * determines the behavior of the interpolation return value. The default is
     464             :  * just interpolation, if `Tag` has prefix  `::Tags::Multiplies` or `Tags::Du`
     465             :  * (this case), the interpolator performs the additional multiplication or time
     466             :  * derivative as a step in the interpolation procedure.
     467             :  */
     468             : template <typename VectorTypeToInterpolate, typename Tag>
     469           1 : struct ScriPlusInterpolationManager<VectorTypeToInterpolate, Tags::Du<Tag>> {
     470             :  public:
     471           0 :   ScriPlusInterpolationManager() = default;
     472           0 :   ScriPlusInterpolationManager(
     473             :       const size_t target_number_of_points, const size_t vector_size,
     474             :       std::unique_ptr<intrp::SpanInterpolator> interpolator)
     475             :       : argument_interpolation_manager_{target_number_of_points, vector_size,
     476             :                                         std::move(interpolator)} {}
     477             : 
     478             :   /// \brief provide data to the interpolation manager.
     479             :   ///
     480             :   /// \details `u_bondi` is a vector of
     481             :   /// inertial times, and `to_interpolate_argument` is the vector of values that
     482             :   /// will be interpolated to target times and differentiated.
     483           1 :   void insert_data(const DataVector& u_bondi,
     484             :                    const VectorTypeToInterpolate& to_interpolate_argument) {
     485             :     argument_interpolation_manager_.insert_data(u_bondi,
     486             :                                                 to_interpolate_argument);
     487             :   }
     488             : 
     489             :   /// \brief Request a target time to be interpolated to when enough data has
     490             :   /// been accumulated.
     491             :   ///
     492             :   /// For optimization, we assume that these are inserted in ascending order.
     493           1 :   void insert_target_time(const double time) {
     494             :     argument_interpolation_manager_.insert_target_time(time);
     495             :   }
     496             : 
     497             :   /// \brief Determines whether enough data before and after the first time in
     498             :   /// the target time queue has been provided to interpolate.
     499             :   ///
     500             :   /// \details If possible, this function will require that the target time to
     501             :   /// be interpolated is reasonably centered on the range, but will settle for
     502             :   /// non-centered data if the time is too early for the given data, which is
     503             :   /// necessary to get the first couple of times out of the simulation. This
     504             :   /// function always returns false if all of the provided data is earlier than
     505             :   /// the next target time, indicating that the caller should wait until more
     506             :   /// data has been provided before interpolating.
     507           1 :   bool first_time_is_ready_to_interpolate() const {
     508             :     return argument_interpolation_manager_.first_time_is_ready_to_interpolate();
     509             :   }
     510             : 
     511           0 :   const std::deque<std::pair<double, double>>& get_u_bondi_ranges() const {
     512             :     return argument_interpolation_manager_.get_u_bondi_ranges();
     513             :   }
     514             : 
     515             :   /// \brief Interpolate to the first target time in the queue, returning both
     516             :   /// the time and the interpolated data at that time.
     517             :   ///
     518             :   /// \note If this function is not able to interpolate to the full accuracy
     519             :   /// using a centered stencil, it will perform the best interpolation
     520             :   /// available.
     521             :   ///
     522             :   /// \warning This function will extrapolate if the target times are
     523             :   /// outside the range of data the class has been provided. This is intentional
     524             :   /// to support small extrapolation at the end of a simulation when no further
     525             :   /// data is available, but for full accuracy, check
     526             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     527             :   /// functions
     528           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_first_time();
     529             : 
     530             :   /// \brief Interpolate to the first target time in the queue, returning both
     531             :   /// the time and the interpolated data at that time, and remove the first time
     532             :   /// from the queue.
     533             :   ///
     534             :   /// \note If this function is not able to interpolate to the full accuracy
     535             :   /// using a centered stencil, it will perform the best interpolation
     536             :   /// available.
     537             :   ///
     538             :   /// \warning This function will extrapolate if the target times are
     539             :   /// outside the range of data the class has been provided. This is intentional
     540             :   /// to support small extrapolation at the end of a simulation when no further
     541             :   /// data is available, but for full accuracy, check
     542             :   /// `first_time_is_ready_to_interpolate` before calling the interpolation
     543             :   /// functions
     544           1 :   std::pair<double, VectorTypeToInterpolate> interpolate_and_pop_first_time();
     545             : 
     546             :   /// \brief return the number of times in the target times queue
     547           1 :   size_t number_of_target_times() const {
     548             :     return argument_interpolation_manager_.number_of_target_times();
     549             :   }
     550             : 
     551             :   /// \brief return the number of data points that have been provided to the
     552             :   /// interpolation manager
     553           1 :   size_t number_of_data_points() const {
     554             :     return argument_interpolation_manager_.number_of_data_points();
     555             :   }
     556             : 
     557             :   // NOLINTNEXTLINE(google-runtime-references)
     558           0 :   void pup(PUP::er& p) { p | argument_interpolation_manager_; }
     559             : 
     560             :  private:
     561             :   // to avoid code duplication, most of the details are stored in an
     562             :   // interpolation manager for the argument tag, and this class is a friend of
     563             :   // the argument tag class
     564             :   /// \cond
     565             :   ScriPlusInterpolationManager<VectorTypeToInterpolate, Tag>
     566             :       argument_interpolation_manager_;
     567             :   /// \endcond
     568             : };
     569             : 
     570             : template <typename VectorTypeToInterpolate, typename Tag>
     571             : std::pair<double, VectorTypeToInterpolate> ScriPlusInterpolationManager<
     572             :     VectorTypeToInterpolate, Tags::Du<Tag>>::interpolate_first_time() {
     573             :   const size_t target_number_of_points =
     574             :       argument_interpolation_manager_.target_number_of_points_;
     575             :   if (argument_interpolation_manager_.target_times_.empty()) {
     576             :     ERROR("There are no target times to interpolate.");
     577             :   }
     578             :   if (argument_interpolation_manager_.to_interpolate_values_.size() <
     579             :       2 * argument_interpolation_manager_.target_number_of_points_) {
     580             :     ERROR("Insufficient data points to continue interpolation: have "
     581             :           << argument_interpolation_manager_.to_interpolate_values_.size()
     582             :           << ", need at least"
     583             :           << 2 * argument_interpolation_manager_.target_number_of_points_);
     584             :   }
     585             :   // note that because we demand at least a certain number before and at least
     586             :   // a certain number after, we are likely to have a surfeit of points for the
     587             :   // interpolator, but this should not cause significant trouble for a
     588             :   // reasonable method.
     589             :   VectorTypeToInterpolate result{argument_interpolation_manager_.vector_size_};
     590             : 
     591             :   VectorTypeToInterpolate interpolation_values{2 * target_number_of_points};
     592             :   VectorTypeToInterpolate lobatto_collocation_values{2 *
     593             :                                                      target_number_of_points};
     594             :   VectorTypeToInterpolate derivative_lobatto_collocation_values{
     595             :       2 * target_number_of_points};
     596             :   DataVector interpolation_times{2 * target_number_of_points};
     597             :   DataVector collocation_points =
     598             :       Spectral::collocation_points<Spectral::Basis::Legendre,
     599             :                                    Spectral::Quadrature::GaussLobatto>(
     600             :           2 * target_number_of_points);
     601             : 
     602             :   const size_t interpolation_data_size =
     603             :       argument_interpolation_manager_.to_interpolate_values_.size();
     604             : 
     605             :   for (size_t i = 0; i < argument_interpolation_manager_.vector_size_; ++i) {
     606             :     // binary search assumes times placed in sorted order
     607             :     auto upper_bound_offset = static_cast<size_t>(std::distance(
     608             :         argument_interpolation_manager_.u_bondi_values_.begin(),
     609             :         std::upper_bound(
     610             :             argument_interpolation_manager_.u_bondi_values_.begin(),
     611             :             argument_interpolation_manager_.u_bondi_values_.end(),
     612             :             argument_interpolation_manager_.target_times_.front(),
     613             :             [&i](const double rhs, const DataVector& lhs) {
     614             :               return rhs < lhs[i];
     615             :             })));
     616             :     size_t lower_bound_offset =
     617             :         upper_bound_offset == 0 ? 0 : upper_bound_offset - 1;
     618             : 
     619             :     if (upper_bound_offset + target_number_of_points >
     620             :         interpolation_data_size) {
     621             :       upper_bound_offset = interpolation_data_size;
     622             :       lower_bound_offset =
     623             :           interpolation_data_size - 2 * target_number_of_points;
     624             :     } else if (lower_bound_offset < target_number_of_points - 1) {
     625             :       lower_bound_offset = 0;
     626             :       upper_bound_offset = 2 * target_number_of_points;
     627             :     } else {
     628             :       lower_bound_offset = lower_bound_offset + 1 - target_number_of_points;
     629             :       upper_bound_offset = lower_bound_offset + 2 * target_number_of_points;
     630             :     }
     631             :     auto interpolation_values_begin =
     632             :         argument_interpolation_manager_.to_interpolate_values_.begin() +
     633             :         static_cast<ptrdiff_t>(lower_bound_offset);
     634             :     auto interpolation_times_begin =
     635             :         argument_interpolation_manager_.u_bondi_values_.begin() +
     636             :         static_cast<ptrdiff_t>(lower_bound_offset);
     637             :     auto interpolation_times_end =
     638             :         argument_interpolation_manager_.u_bondi_values_.begin() +
     639             :         static_cast<ptrdiff_t>(upper_bound_offset);
     640             : 
     641             :     // interpolate using the data sets in the restricted iterators
     642             :     auto value_it = interpolation_values_begin;
     643             :     size_t vector_position = 0;
     644             :     for (auto time_it = interpolation_times_begin;
     645             :          time_it != interpolation_times_end;
     646             :          ++time_it, ++value_it, ++vector_position) {
     647             :       interpolation_values[vector_position] = (*value_it)[i];
     648             :       interpolation_times[vector_position] = (*time_it)[i];
     649             :     }
     650             :     for (size_t j = 0; j < lobatto_collocation_values.size(); ++j) {
     651             :       lobatto_collocation_values[j] =
     652             :           argument_interpolation_manager_.interpolator_->interpolate(
     653             :               gsl::span<const double>(interpolation_times.data(),
     654             :                                       interpolation_times.size()),
     655             :               gsl::span<const typename VectorTypeToInterpolate::value_type>(
     656             :                   interpolation_values.data(), interpolation_values.size()),
     657             :               // affine transformation between the Gauss-Lobatto collocation
     658             :               // points and the physical times
     659             :               (collocation_points[j] + 1.0) * 0.5 *
     660             :                       (interpolation_times[interpolation_times.size() - 1] -
     661             :                        interpolation_times[0]) +
     662             :                   interpolation_times[0]);
     663             :     }
     664             :     // note the coordinate transformation to and from the Gauss-Lobatto basis
     665             :     // range [-1, 1]
     666             :     apply_matrices(
     667             :         make_not_null(&derivative_lobatto_collocation_values),
     668             :         make_array<1>(
     669             :             Spectral::differentiation_matrix<
     670             :                 Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>(
     671             :                 lobatto_collocation_values.size())),
     672             :         lobatto_collocation_values,
     673             :         Index<1>(lobatto_collocation_values.size()));
     674             : 
     675             :     result[i] =
     676             :         argument_interpolation_manager_.interpolator_->interpolate(
     677             :             gsl::span<const double>(collocation_points.data(),
     678             :                                     collocation_points.size()),
     679             :             gsl::span<const typename VectorTypeToInterpolate::value_type>(
     680             :                 derivative_lobatto_collocation_values.data(),
     681             :                 derivative_lobatto_collocation_values.size()),
     682             :             2.0 *
     683             :                     (argument_interpolation_manager_.target_times_.front() -
     684             :                      interpolation_times[0]) /
     685             :                     (interpolation_times[interpolation_times.size() - 1] -
     686             :                      interpolation_times[0]) -
     687             :                 1.0) *
     688             :         2.0 /
     689             :         (interpolation_times[interpolation_times.size() - 1] -
     690             :          interpolation_times[0]);
     691             :   }
     692             :   return std::make_pair(argument_interpolation_manager_.target_times_.front(),
     693             :                         std::move(result));
     694             : }
     695             : 
     696             : template <typename VectorTypeToInterpolate, typename Tag>
     697             : std::pair<double, VectorTypeToInterpolate> ScriPlusInterpolationManager<
     698             :     VectorTypeToInterpolate, Tags::Du<Tag>>::interpolate_and_pop_first_time() {
     699             :   std::pair<double, VectorTypeToInterpolate> interpolated =
     700             :       interpolate_first_time();
     701             :   argument_interpolation_manager_.target_times_.pop_front();
     702             : 
     703             :   if (not argument_interpolation_manager_.target_times_.empty()) {
     704             :     argument_interpolation_manager_.remove_unneeded_early_times();
     705             :   }
     706             :   return interpolated;
     707             : }
     708             : }  // namespace Cce

Generated by: LCOV version 1.14