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