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