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