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