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