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