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