ReadBoundaryDataH5.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <boost/iterator/zip_iterator.hpp>
8 #include <boost/tuple/tuple.hpp>
9 #include <boost/tuple/tuple_comparison.hpp>
10 #include <memory>
11 #include <string>
12 
16 #include "ErrorHandling/Assert.hpp"
17 #include "Evolution/Systems/Cce/BoundaryData.hpp"
18 #include "Evolution/Systems/Cce/SpecBoundaryData.hpp"
19 #include "IO/H5/Dat.hpp"
20 #include "IO/H5/File.hpp"
21 #include "IO/H5/Version.hpp"
22 #include "NumericalAlgorithms/Interpolation/SpanInterpolator.hpp"
24 
25 namespace Cce {
26 namespace Tags {
27 namespace detail {
28 // tags for use in the buffers for the modal input worldtube data management
29 // classes
30 using SpatialMetric =
34 
35 // radial derivative prefix tag to be used with the modal input worldtube data
36 template <typename Tag>
37 struct Dr : db::SimpleTag, db::PrefixTag {
38  using type = typename Tag::type;
39  using tag = Tag;
40  static std::string name() noexcept { return "Dr(" + Tag::name() + ")"; }
41 };
42 
43 // tag for the string for accessing the quantity associated with `Tag` in
44 // worldtube h5 file
45 template <typename Tag>
46 struct InputDataSet : db::SimpleTag, db::PrefixTag {
47  using type = std::string;
48  using tag = Tag;
49  static std::string name() noexcept {
50  return "InputDataSet(" + Tag::name() + ")";
51  }
52 };
53 } // namespace detail
54 } // namespace Tags
55 
56 namespace detail {
57 // generates the component dataset name in the worldtube file based on the
58 // tensor indices requested. For instance, if called with arguments ("/g", 0,1),
59 // it returns the dataset name "/gxy".
60 template <typename... T>
61 std::string dataset_name_for_component(std::string base_name,
62  const T... indices) noexcept {
63  const auto add_index = [&base_name](size_t index) noexcept {
64  ASSERT(index < 3, "The character-arithmetic index must be less than 3.");
65  base_name += static_cast<char>('x' + index);
66  };
67  EXPAND_PACK_LEFT_TO_RIGHT(add_index(indices));
68  // void cast so that compilers can tell it's used.
69  (void)add_index;
70  return base_name;
71 }
72 
73 // the full set of tensors to be extracted from the worldtube h5 file
74 using cce_input_tags = tmpl::list<
75  Tags::detail::SpatialMetric, Tags::detail::Dr<Tags::detail::SpatialMetric>,
76  ::Tags::dt<Tags::detail::SpatialMetric>, Tags::detail::Shift,
77  Tags::detail::Dr<Tags::detail::Shift>, ::Tags::dt<Tags::detail::Shift>,
78  Tags::detail::Lapse, Tags::detail::Dr<Tags::detail::Lapse>,
80 
81 using reduced_cce_input_tags =
82  tmpl::list<Spectral::Swsh::Tags::SwshTransform<Tags::BondiBeta>,
91 
92 // creates a pair of indices such that the difference is `2 *
93 // interpolator_length + pad`, centered around `time`, and bounded by
94 // `lower_bound` and `upper_bound`. If it cannot be centered, it gives a span
95 // that is appropriately sized and bounded by the supplied bounds. If the bounds
96 // are too constraining for the necessary size, it gives a span that is the
97 // correct size starting at `lower bound`, but not constrained by `upper_bound`
98 std::pair<size_t, size_t> create_span_for_time_value(
99  double time, size_t pad, size_t interpolator_length, size_t lower_bound,
100  size_t upper_bound, const DataVector& time_buffer) noexcept;
101 } // namespace detail
102 
103 /// \cond
104 class SpecWorldtubeH5BufferUpdater;
105 /// \endcond
106 
107 /*!
108  * \brief Abstract base class for utilities that are able to perform the buffer
109  * updating procedure needed by the `WorldtubeDataManager`.
110  *
111  * \details The methods that are required to be overridden in the derived
112  * classes are:
113  * - `WorldtubeBufferUpdater::update_buffers_for_time()`:
114  * updates the buffers passed by pointer and the `time_span_start` and
115  * `time_span_end` to be appropriate for the requested `time`,
116  * `interpolator_length`, and `buffer_depth`.
117  * - `WorldtubeBufferUpdater::get_clone()`
118  * clone function to obtain a `std::unique_ptr` of the base
119  * `WorldtubeBufferUpdater`, needed to pass around the factory-created
120  * object.
121  * - `WorldtubeBufferUpdater::time_is_outside_range()`
122  * the override should return `true` if the `time` could be used in a
123  * `update_buffers_for_time` call given the data available to the derived
124  * class, and `false` otherwise
125  * - `WorldtubeBufferUpdater::get_l_max()`
126  * The override should return the `l_max` it uses in the
127  * Goldberg modal data placed in the buffers.
128  * - `WorldtubeBufferUpdater::get_extraction_radius()`
129  * The override should return the coordinate radius associated with the modal
130  * worldtube data that it supplies in the buffer update function. This is
131  * currently assumed to be a single double, but may be generalized in future
132  * to be time-dependent.
133  * - `WorldtubeBufferUpdater::get_time_buffer`
134  * The override should return the vector of times that it can produce modal
135  * data at. For instance, if associated with a file input, this will be the
136  * times at each of the rows of the time-series data.
137  */
138 class WorldtubeBufferUpdater : public PUP::able {
139  using creatable_classes = tmpl::list<SpecWorldtubeH5BufferUpdater>;
140 
142 
143  virtual double update_buffers_for_time(
144  gsl::not_null<Variables<detail::cce_input_tags>*> buffers,
145  gsl::not_null<size_t*> time_span_start,
146  gsl::not_null<size_t*> time_span_end, double time,
147  size_t interpolator_length, size_t buffer_depth) const noexcept = 0;
148 
149  virtual std::unique_ptr<WorldtubeBufferUpdater> get_clone() const
150  noexcept = 0;
151 
152  virtual bool time_is_outside_range(double time) const noexcept = 0;
153 
154  virtual size_t get_l_max() const noexcept = 0;
155 
156  virtual double get_extraction_radius() const noexcept = 0;
157 
158  virtual bool radial_derivatives_need_renormalization() const noexcept = 0;
159 
160  virtual DataVector& get_time_buffer() noexcept = 0;
161 };
162 
163 /// A `WorldtubeBufferUpdater` specialized to the CCE input worldtube H5 file
164 /// produced by SpEC.
166  public:
167  // charm needs the empty constructor
168  SpecWorldtubeH5BufferUpdater() = default;
169 
170  /// The constructor takes the filename of the SpEC h5 file that will be used
171  /// for boundary data. Note that this assumes that the input data has
172  /// correctly-normalized radial derivatives, and that the extraction radius is
173  /// encoded as an integer in the filename.
175  const std::string& cce_data_filename) noexcept;
176 
178 
179  explicit SpecWorldtubeH5BufferUpdater(CkMigrateMessage* /*unused*/) noexcept {
180  }
181 
182  /// update the `buffers`, `time_span_start`, and `time_span_end` with
183  /// time-varies-fastest, Goldberg modal data and the start and end index in
184  /// the member `time_buffer_` covered by the newly updated `buffers`. The
185  /// function returns the next time at which a full update will occur. If
186  /// called again at times earlier than the next full update time, it will
187  /// leave the `buffers` unchanged and again return the next needed time.
188  double update_buffers_for_time(
189  gsl::not_null<Variables<detail::cce_input_tags>*> buffers,
190  gsl::not_null<size_t*> time_span_start,
191  gsl::not_null<size_t*> time_span_end, double time,
192  size_t interpolator_length, size_t buffer_depth) const noexcept override;
193 
194  std::unique_ptr<WorldtubeBufferUpdater> get_clone() const noexcept override;
195 
196  /// The time can only be supported in the buffer update if it is between the
197  /// first and last time of the input file.
198  bool time_is_outside_range(double time) const noexcept override;
199 
200  /// retrieves the l_max of the input file
201  size_t get_l_max() const noexcept override { return l_max_; }
202 
203  /// retrieves the extraction radius encoded in the filename
204  double get_extraction_radius() const noexcept override {
205  return extraction_radius_;
206  }
207 
208  /// The time buffer is supplied by non-const reference to allow views to
209  /// easily point into the buffer.
210  ///
211  /// \warning Altering this buffer outside of the constructor of this class
212  /// results in undefined behavior! This should be supplied by const reference
213  /// once there is a convenient method of producing a const view of a vector
214  /// type.
215  DataVector& get_time_buffer() noexcept override { return time_buffer_; }
216 
217  bool radial_derivatives_need_renormalization() const noexcept override {
218  return radial_derivatives_need_renormalization_;
219  }
220 
221  /// Serialization for Charm++.
222  void pup(PUP::er& p) noexcept override;
223 
224  private:
225  void update_buffer(gsl::not_null<ComplexModalVector*> buffer_to_update,
226  const h5::Dat& read_data, size_t time_span_start,
227  size_t time_span_end) const noexcept;
228 
229  bool radial_derivatives_need_renormalization_ = false;
230  double extraction_radius_ = 1.0;
231  size_t l_max_ = 0;
232 
234  std::string filename_;
235 
236  tuples::tagged_tuple_from_typelist<
238  dataset_names_;
239 
240  // stores all the times in the input file
241  DataVector time_buffer_;
242 };
243 
244 /// \cond
246 /// \endcond
247 
248 class ReducedWorldtubeBufferUpdater : public PUP::able {
249  using creatable_classes = tmpl::list<ReducedSpecWorldtubeH5BufferUpdater>;
250 
252 
253  virtual double update_buffers_for_time(
254  gsl::not_null<Variables<detail::reduced_cce_input_tags>*> buffers,
255  gsl::not_null<size_t*> time_span_start,
256  gsl::not_null<size_t*> time_span_end, double time,
257  size_t interpolator_length, size_t buffer_depth) const noexcept = 0;
258 
259  virtual std::unique_ptr<ReducedWorldtubeBufferUpdater> get_clone() const
260  noexcept = 0;
261 
262  virtual bool time_is_outside_range(double time) const noexcept = 0;
263 
264  virtual size_t get_l_max() const noexcept = 0;
265 
266  virtual double get_extraction_radius() const noexcept = 0;
267 
268  virtual DataVector& get_time_buffer() noexcept = 0;
269 };
270 
271 /// A `WorldtubeBufferUpdater` specialized to the CCE input worldtube H5 file
272 /// produced by the reduced SpEC format.
275  public:
276  // charm needs the empty constructor
278 
279  /// The constructor takes the filename of the SpEC h5 file that will be used
280  /// for boundary data. Note that this assumes that the input data has
281  /// correctly-normalized radial derivatives, and that the extraction radius is
282  /// encoded as an integer in the filename.
284  const std::string& cce_data_filename) noexcept;
285 
287 
289  CkMigrateMessage* /*unused*/) noexcept {}
290 
291  /// update the `buffers`, `time_span_start`, and `time_span_end` with
292  /// time-varies-fastest, Goldberg modal data and the start and end index in
293  /// the member `time_buffer_` covered by the newly updated `buffers`.
294  double update_buffers_for_time(
295  gsl::not_null<Variables<detail::reduced_cce_input_tags>*> buffers,
296  gsl::not_null<size_t*> time_span_start,
297  gsl::not_null<size_t*> time_span_end, double time,
298  size_t interpolator_length, size_t buffer_depth) const noexcept override;
299 
301  noexcept override {
302  return std::make_unique<ReducedSpecWorldtubeH5BufferUpdater>(filename_);
303  }
304 
305  /// The time can only be supported in the buffer update if it is between the
306  /// first and last time of the input file.
307  bool time_is_outside_range(const double time) const noexcept override {
308  return time < time_buffer_[0] or
309  time > time_buffer_[time_buffer_.size() - 1];
310  }
311 
312  /// retrieves the l_max of the input file
313  size_t get_l_max() const noexcept override { return l_max_; }
314 
315  /// retrieves the extraction radius encoded in the filename
316  double get_extraction_radius() const noexcept override {
317  return extraction_radius_;
318  }
319 
320  /// The time buffer is supplied by non-const reference to allow views to
321  /// easily point into the buffer.
322  ///
323  /// \warning Altering this buffer outside of the constructor of this class
324  /// results in undefined behavior! This should be supplied by const reference
325  /// once there is a convenient method of producing a const view of a vector
326  /// type.
327  DataVector& get_time_buffer() noexcept override { return time_buffer_; }
328 
329  /// Serialization for Charm++.
330  void pup(PUP::er& p) noexcept override;
331 
332  private:
333  void update_buffer(gsl::not_null<ComplexModalVector*> buffer_to_update,
334  const h5::Dat& read_data, const size_t& time_span_start,
335  const size_t& time_span_end) const noexcept;
336 
337  double extraction_radius_ = 1.0;
338  size_t l_max_ = 0;
339 
341  std::string filename_;
342 
343  tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
344  Tags::detail::InputDataSet, detail::reduced_cce_input_tags>>
345  dataset_names_;
346 
347  // stores all the times in the input file
348  DataVector time_buffer_;
349 };
350 
351 /*!
352  * \brief Manages the cached buffer data associated with a CCE worldtube and
353  * interpolates to requested time points to provide worldtube boundary data to
354  * the main evolution routines.
355  *
356  * \details The maintained buffer will be maintained at a length that is set by
357  * the `Interpolator` and the `buffer_depth` also passed to the constructor. A
358  * longer depth will ensure that the buffer updater is called less frequently,
359  * which is useful for slow updaters (e.g. those that perform file access).
360  * The main functionality is provided by the
361  * `WorldtubeDataManager::populate_hypersurface_boundary_data()` member
362  * function that handles buffer updating and boundary computation.
363  */
365  public:
366  // charm needs an empty constructor.
367  WorldtubeDataManager() noexcept = default;
368 
371  const size_t l_max, const size_t buffer_depth,
372  std::unique_ptr<intrp::SpanInterpolator> interpolator) noexcept
373  : buffer_updater_{std::move(buffer_updater)},
374  l_max_{l_max},
375  interpolated_coefficients_{
377  buffer_depth_{buffer_depth},
378  interpolator_{std::move(interpolator)} {
379  if (UNLIKELY(
380  buffer_updater_->get_time_buffer().size() <
381  2 * interpolator_->required_number_of_points_before_and_after() +
382  buffer_depth)) {
383  ERROR(
384  "The specified buffer updater doesn't have enough time points to "
385  "supply the requested interpolation buffer. This almost certainly "
386  "indicates that the corresponding file hasn't been created properly, "
387  "but might indicate that the `buffer_depth` template parameter is "
388  "too large or the specified Interpolator requests too many points");
389  }
390 
391  const size_t size_of_buffer =
392  square(l_max + 1) *
393  (buffer_depth +
394  2 * interpolator_->required_number_of_points_before_and_after());
395  coefficients_buffers_ = Variables<detail::cce_input_tags>{size_of_buffer};
396  }
397 
398  /*!
399  * \brief Update the `boundary_data_box` entries for all tags in
400  * `Tags::characteristic_worldtube_boundary_tags` to the boundary data at
401  * `time`.
402  *
403  * \details First, if the stored buffer requires updating, it will be updated
404  * via the `buffer_updater_` supplied in the constructor. Then, each of the
405  * spatial metric, shift, lapse, and each of their radial and time derivatives
406  * are interpolated across buffer points to the requested time value (via the
407  * `Interpolator` provided in the constructor). Finally, that data is supplied
408  * to the `create_bondi_boundary_data()`, which updates the
409  * `boundary_data_box` with the Bondi spin-weighted scalars determined from
410  * the interpolated Cartesian data.
411  *
412  * Returns `true` if the time can be supplied from the `buffer_updater_`, and
413  * `false` otherwise. No tags are updated if `false` is returned.
414  */
415  template <typename TagList>
417  const gsl::not_null<db::DataBox<TagList>*> boundary_data_box,
418  const double time) noexcept {
419  if (buffer_updater_->time_is_outside_range(time)) {
420  return false;
421  }
422  buffer_updater_->update_buffers_for_time(
423  make_not_null(&coefficients_buffers_), make_not_null(&time_span_start_),
424  make_not_null(&time_span_end_), time,
425  interpolator_->required_number_of_points_before_and_after(),
426  buffer_depth_);
427  const auto interpolation_time_span = detail::create_span_for_time_value(
428  time, 0, interpolator_->required_number_of_points_before_and_after(),
429  time_span_start_, time_span_end_, buffer_updater_->get_time_buffer());
430 
431  // search through and find the two interpolation points the time point is
432  // between. If we can, put the range for the interpolation centered on the
433  // desired point. If that can't be done (near the start or the end of the
434  // simulation), make the range terminated at the start or end of the cached
435  // data and extending for the desired range in the other direction.
436  const size_t buffer_span_size = time_span_end_ - time_span_start_;
437  const size_t interpolation_span_size =
438  interpolation_time_span.second - interpolation_time_span.first;
439 
440  const DataVector time_points{buffer_updater_->get_time_buffer().data() +
441  interpolation_time_span.first,
442  interpolation_span_size};
443 
444  auto interpolate_from_column = [
445  &time, &time_points, &buffer_span_size, &interpolation_time_span,
446  &interpolation_span_size,
447  this
448  ](auto data, const size_t column) noexcept {
449  auto interp_val = interpolator_->interpolate(
450  gsl::span<const double>(time_points.data(), time_points.size()),
452  data + column * buffer_span_size +
453  (interpolation_time_span.first - time_span_start_),
454  interpolation_span_size),
455  time);
456  return interp_val;
457  };
458 
459  // the ComplexModalVectors should be provided from the buffer_updater_ in
460  // 'Goldberg' format, so we iterate over modes and convert to libsharp
461  // format.
462 
463  // we'll just use this buffer to reference into the actual data to satisfy
464  // the swsh interface requirement that the spin-weight be labelled with
465  // `SpinWeighted`
466  SpinWeighted<ComplexModalVector, 0> spin_weighted_buffer;
467  for (const auto& libsharp_mode :
469  for (size_t i = 0; i < 3; ++i) {
470  for (size_t j = i; j < 3; ++j) {
471  tmpl::for_each<
472  tmpl::list<Tags::detail::SpatialMetric,
473  Tags::detail::Dr<Tags::detail::SpatialMetric>,
475  this, &i, &j, &libsharp_mode, &interpolate_from_column, &
476  spin_weighted_buffer
477  ](auto tag_v) noexcept {
478  using tag = typename decltype(tag_v)::type;
479  spin_weighted_buffer.set_data_ref(
480  get<tag>(interpolated_coefficients_).get(i, j).data(),
482  Spectral::Swsh::goldberg_modes_to_libsharp_modes_single_pair(
483  libsharp_mode, make_not_null(&spin_weighted_buffer), 0,
484  interpolate_from_column(
485  get<tag>(coefficients_buffers_).get(i, j).data(),
487  l_max_, libsharp_mode.l,
488  static_cast<int>(libsharp_mode.m))),
489  interpolate_from_column(
490  get<tag>(coefficients_buffers_).get(i, j).data(),
492  l_max_, libsharp_mode.l,
493  -static_cast<int>(libsharp_mode.m))));
494  });
495  }
496  tmpl::for_each<tmpl::list<Tags::detail::Shift,
497  Tags::detail::Dr<Tags::detail::Shift>,
499  this, &i, &libsharp_mode, &interpolate_from_column, &
500  spin_weighted_buffer
501  ](auto tag_v) noexcept {
502  using tag = typename decltype(tag_v)::type;
503  spin_weighted_buffer.set_data_ref(
504  get<tag>(interpolated_coefficients_).get(i).data(),
506  Spectral::Swsh::goldberg_modes_to_libsharp_modes_single_pair(
507  libsharp_mode, make_not_null(&spin_weighted_buffer), 0,
508  interpolate_from_column(
509  get<tag>(coefficients_buffers_).get(i).data(),
511  l_max_, libsharp_mode.l,
512  static_cast<int>(libsharp_mode.m))),
513  interpolate_from_column(
514  get<tag>(coefficients_buffers_).get(i).data(),
516  l_max_, libsharp_mode.l,
517  -static_cast<int>(libsharp_mode.m))));
518  });
519  }
520  tmpl::for_each<
521  tmpl::list<Tags::detail::Lapse, Tags::detail::Dr<Tags::detail::Lapse>,
523  this, &libsharp_mode, &interpolate_from_column, &spin_weighted_buffer
524  ](auto tag_v) noexcept {
525  using tag = typename decltype(tag_v)::type;
526  spin_weighted_buffer.set_data_ref(
527  get(get<tag>(interpolated_coefficients_)).data(),
529  Spectral::Swsh::goldberg_modes_to_libsharp_modes_single_pair(
530  libsharp_mode, make_not_null(&spin_weighted_buffer), 0,
531  interpolate_from_column(get(get<tag>(coefficients_buffers_)).data(),
533  l_max_, libsharp_mode.l,
534  static_cast<int>(libsharp_mode.m))),
535  interpolate_from_column(get(get<tag>(coefficients_buffers_)).data(),
537  l_max_, libsharp_mode.l,
538  -static_cast<int>(libsharp_mode.m))));
539  });
540  }
541  // At this point, we have a collection of 9 tensors of libsharp
542  // coefficients. This is what the boundary data calculation utility takes
543  // as an input, so we now hand off the control flow to the boundary and
544  // gauge transform utility
545  if (buffer_updater_->radial_derivatives_need_renormalization()) {
547  boundary_data_box,
548  get<Tags::detail::SpatialMetric>(interpolated_coefficients_),
550  interpolated_coefficients_),
551  get<Tags::detail::Dr<Tags::detail::SpatialMetric>>(
552  interpolated_coefficients_),
553  get<Tags::detail::Shift>(interpolated_coefficients_),
555  interpolated_coefficients_),
556  get<Tags::detail::Dr<Tags::detail::Shift>>(
557  interpolated_coefficients_),
558  get<Tags::detail::Lapse>(interpolated_coefficients_),
560  interpolated_coefficients_),
561  get<Tags::detail::Dr<Tags::detail::Lapse>>(
562  interpolated_coefficients_),
563  buffer_updater_->get_extraction_radius(), l_max_);
564  } else {
566  boundary_data_box,
567  get<Tags::detail::SpatialMetric>(interpolated_coefficients_),
569  interpolated_coefficients_),
570  get<Tags::detail::Dr<Tags::detail::SpatialMetric>>(
571  interpolated_coefficients_),
572  get<Tags::detail::Shift>(interpolated_coefficients_),
574  interpolated_coefficients_),
575  get<Tags::detail::Dr<Tags::detail::Shift>>(
576  interpolated_coefficients_),
577  get<Tags::detail::Lapse>(interpolated_coefficients_),
579  interpolated_coefficients_),
580  get<Tags::detail::Dr<Tags::detail::Lapse>>(
581  interpolated_coefficients_),
582  buffer_updater_->get_extraction_radius(), l_max_);
583  }
584  return true;
585  }
586 
587  /// retrieves the l_max that will be supplied to the \ref DataBoxGroup in
588  /// `populate_hypersurface_boundary_data()`
589  size_t get_l_max() const noexcept { return l_max_; }
590 
591  /// retrieves the current time span associated with the `buffer_updater_` for
592  /// diagnostics
594  return std::make_pair(time_span_start_, time_span_end_);
595  }
596 
597  /// Serialization for Charm++.
598  void pup(PUP::er& p) noexcept { // NOLINT
599  p | time_span_start_;
600  p | time_span_end_;
601  p | l_max_;
602  p | buffer_depth_;
603  p | interpolator_;
604  if (p.isUnpacking()) {
605  const size_t size_of_buffer =
606  square(l_max_ + 1) *
607  (buffer_depth_ +
608  2 * interpolator_->required_number_of_points_before_and_after());
609  coefficients_buffers_ = Variables<detail::cce_input_tags>{size_of_buffer};
610  interpolated_coefficients_ = Variables<detail::cce_input_tags>{
612  }
613  }
614 
615  private:
617  size_t time_span_start_ = 0;
618  size_t time_span_end_ = 0;
619  size_t l_max_ = 0;
620 
621  // These buffers are just kept around to avoid allocations; they're
622  // updated every time a time is requested
623  Variables<detail::cce_input_tags> interpolated_coefficients_;
624 
625  // note: buffers store data in a 'time-varies-fastest' manner
626  Variables<detail::cce_input_tags> coefficients_buffers_;
627 
628  size_t buffer_depth_ = 0;
629 
631 };
632 
633 /*!
634  * \brief Manages the 'reduced' cached buffer dataset associated with a CCE
635  * worldtube and interpolates to requested time points to provide worldtube
636  * boundary data to the main evolution routines.
637  *
638  * \details The maintained buffer will be kept at a length that is set by
639  * the `Interpolator` and the `buffer_depth` also passed to the constructor. A
640  * longer depth will ensure that the buffer updater is called less frequently,
641  * which is useful for slow updaters (e.g. those that perform file access).
642  * The main functionality is provided by the
643  * `WorldtubeDataManager::populate_hypersurface_boundary_data()` member
644  * function that handles buffer updating and boundary computation. This version
645  * of the data manager handles the 9 scalars of
646  * `detail::reduced_cce_input_tags`, rather than direct metric components
647  * handled by `WorldtubeDataManager`. The set of 9 scalars is a far leaner
648  * (factor of ~4) data storage format.
649  */
651  public:
652  // charm needs an empty constructor.
653  ReducedWorldtubeDataManager() = default;
654 
657  size_t l_max, size_t buffer_depth,
658  std::unique_ptr<intrp::SpanInterpolator> interpolator) noexcept;
659 
660  /*!
661  * \brief Update the `boundary_data_box` entries for all tags in
662  * `Tags::characteristic_worldtube_boundary_tags` to the boundary data at
663  * `time`.
664  *
665  * \details First, if the stored buffer requires updating, it will be updated
666  * via the `buffer_updater_` supplied in the constructor. Then, each of the
667  * 9 spin-weighted scalars in `detail::reduced_cce_input_tags`
668  * are interpolated across buffer points to the requested time value (via the
669  * `Interpolator` provided in the constructor). Finally, the remaining two
670  * scalars not directly supplied in the input file are calculated in-line and
671  * put in the \ref DataBoxGroup.
672  *
673  * Returns `true` if the time can be supplied from the `buffer_updater_`, and
674  * `false` otherwise. No tags are updated if `false` is returned.
675  */
676  template <typename TagList>
677  bool populate_hypersurface_boundary_data(
678  gsl::not_null<db::DataBox<TagList>*> boundary_data_box,
679  double time) noexcept;
680 
681  /// retrieves the l_max that will be supplied to the \ref DataBoxGroup in
682  /// `populate_hypersurface_boundary_data()`
683  size_t get_l_max() const noexcept { return l_max_; }
684 
685  /// retrieves the current time span associated with the `buffer_updater_` for
686  /// diagnostics
688  return std::make_pair(time_span_start_, time_span_end_);
689  }
690 
691  /// Serialization for Charm++.
692  void pup(PUP::er& p) noexcept; // NOLINT
693 
694  private:
696  size_t time_span_start_ = 0;
697  size_t time_span_end_ = 0;
698  size_t l_max_ = 0;
699 
700  // These buffers are just kept around to avoid allocations; they're
701  // updated every time a time is requested
702  Variables<detail::reduced_cce_input_tags> interpolated_coefficients_;
703 
704  // note: buffers store data in an 'time-varies-fastest' manner
705  Variables<detail::reduced_cce_input_tags> coefficients_buffers_;
706 
707  size_t buffer_depth_ = 0;
708 
710 };
711 
712 template <typename TagList>
714  const gsl::not_null<db::DataBox<TagList>*> boundary_data_box,
715  const double time) noexcept {
716  if (buffer_updater_->time_is_outside_range(time)) {
717  return false;
718  }
719  buffer_updater_->update_buffers_for_time(
720  make_not_null(&coefficients_buffers_), make_not_null(&time_span_start_),
721  make_not_null(&time_span_end_), time,
722  interpolator_->required_number_of_points_before_and_after(),
723  buffer_depth_);
724  auto interpolation_time_span = detail::create_span_for_time_value(
725  time, 0, interpolator_->required_number_of_points_before_and_after(),
726  time_span_start_, time_span_end_, buffer_updater_->get_time_buffer());
727 
728  // search through and find the two interpolation points the time point is
729  // between. If we can, put the range for the interpolation centered on the
730  // desired point. If that can't be done (near the start or the end of the
731  // simulation), make the range terminated at the start or end of the cached
732  // data and extending for the desired range in the other direction.
733  const size_t buffer_span_size = time_span_end_ - time_span_start_;
734  const size_t interpolation_span_size =
735  interpolation_time_span.second - interpolation_time_span.first;
736 
737  DataVector time_points{
738  buffer_updater_->get_time_buffer().data() + interpolation_time_span.first,
739  interpolation_span_size};
740 
741  auto interpolate_from_column =
742  [&time, &time_points, &buffer_span_size, &interpolation_time_span,
743  &interpolation_span_size, this](auto data, size_t column) {
744  const auto interp_val = interpolator_->interpolate(
745  gsl::span<const double>(time_points.data(), time_points.size()),
747  data + column * (buffer_span_size) +
748  (interpolation_time_span.first - time_span_start_),
749  interpolation_span_size),
750  time);
751  return interp_val;
752  };
753 
754  // the ComplexModalVectors should be provided from the buffer_updater_ in
755  // 'Goldberg' format, so we iterate over modes and convert to libsharp
756  // format.
757 
758  // we'll just use this buffer to reference into the actual data to satisfy
759  // the swsh interface requirement that the spin-weight be labelled with
760  // `SpinWeighted`
761  for (const auto& libsharp_mode :
763  tmpl::for_each<detail::reduced_cce_input_tags>([
764  this, &libsharp_mode, &interpolate_from_column
765  ](auto tag_v) noexcept {
766  using tag = typename decltype(tag_v)::type;
767  Spectral::Swsh::goldberg_modes_to_libsharp_modes_single_pair(
768  libsharp_mode,
769  make_not_null(&get(get<tag>(interpolated_coefficients_))), 0,
770  interpolate_from_column(
771  get(get<tag>(coefficients_buffers_)).data().data(),
773  l_max_, libsharp_mode.l, static_cast<int>(libsharp_mode.m))),
774  interpolate_from_column(
775  get(get<tag>(coefficients_buffers_)).data().data(),
777  l_max_, libsharp_mode.l,
778  -static_cast<int>(libsharp_mode.m))));
779  });
780  }
781 
782  // just inverse transform the 'direct' tags
783  tmpl::for_each<tmpl::transform<detail::reduced_cce_input_tags,
784  tmpl::bind<db::remove_tag_prefix, tmpl::_1>>>(
785  [this, &boundary_data_box](auto tag_v) {
786  using tag = typename decltype(tag_v)::type;
787  db::mutate<Tags::BoundaryValue<tag>>(
788  boundary_data_box,
789  [this](
790  const gsl::not_null<
792  boundary_value,
794  modal_boundary_value) noexcept {
796  l_max_, 1, make_not_null(&get(*boundary_value)),
797  get(modal_boundary_value));
798  },
799  get<Spectral::Swsh::Tags::SwshTransform<tag>>(
800  interpolated_coefficients_));
801  });
802 
803  // there's only a couple of tags desired by the core computation that aren't
804  // stored in the 'reduced' format, so we perform the remaining computation
805  // in-line here.
806 
807  // \partial_u J:
808  db::mutate<Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>>(
809  boundary_data_box,
815  du_r_divided_by_r) noexcept {
816  get(*du_j) = get(h) - get(r) * get(du_r_divided_by_r) * get(dr_j);
817  },
818  db::get<Tags::BoundaryValue<Tags::BondiH>>(*boundary_data_box),
819  db::get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(*boundary_data_box),
820  db::get<Tags::BoundaryValue<Tags::BondiR>>(*boundary_data_box),
821  get<Tags::BoundaryValue<Tags::DuRDividedByR>>(*boundary_data_box));
822 
823  // \partial_r U:
824  db::mutate<Tags::BoundaryValue<Tags::Dr<Tags::BondiU>>>(
825  boundary_data_box,
830  const Scalar<SpinWeighted<ComplexDataVector, 2>>& j) noexcept {
831  // allocation
833  k.data() = sqrt(1.0 + get(j).data() * conj(get(j).data()));
834  get(*dr_u).data() =
835  exp(2.0 * get(beta).data()) / square(get(r).data()) *
836  (k.data() * get(q).data() - get(j).data() * conj(get(q).data()));
837  },
838  get<Tags::BoundaryValue<Tags::BondiBeta>>(*boundary_data_box),
839  get<Tags::BoundaryValue<Tags::BondiR>>(*boundary_data_box),
840  get<Tags::BoundaryValue<Tags::BondiQ>>(*boundary_data_box),
841  get<Tags::BoundaryValue<Tags::BondiJ>>(*boundary_data_box));
842  return true;
843 }
844 } // namespace Cce
constexpr size_t size_of_libsharp_coefficient_vector(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonics coefficients tha...
Definition: SwshCoefficients.hpp:34
Defines class Matrix.
#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: BoundaryData.cpp:21
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:562
Defines class h5::H5File.
#define UNLIKELY(x)
Definition: Gsl.hpp:72
const CoefficientsMetadata & cached_coefficients_metadata(const size_t l_max) noexcept
Generation function for obtaining a CoefficientsMetadata object which is computed by the libsharp cal...
Definition: SwshCoefficients.cpp:32
size_t get_l_max() const noexcept
retrieves the l_max that will be supplied to the DataBox in populate_hypersurface_boundary_data() ...
Definition: ReadBoundaryDataH5.hpp:683
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:64
A WorldtubeBufferUpdater specialized to the CCE input worldtube H5 file produced by the reduced SpEC ...
Definition: ReadBoundaryDataH5.hpp:273
void create_bondi_boundary_data(const gsl::not_null< db::DataBox< DataBoxTagList > *> bondi_boundary_data, const tnsr::iaa< DataVector, 3 > &phi, const tnsr::aa< DataVector, 3 > &pi, const tnsr::aa< DataVector, 3 > &spacetime_metric, const double extraction_radius, const size_t l_max) noexcept
Process the worldtube data from generalized harmonic quantities to desired Bondi quantities, placing the result in the passed DataBox.
Definition: BoundaryData.hpp:784
Defines macros to allow serialization of abstract template base classes.
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new list of Tags by wrapping each tag in TagList using the Wrapper.
Definition: DataBoxTag.hpp:515
Prefix tag representing the spin-weighted spherical harmonic transform of a spin-weighted scalar...
Definition: SwshTags.hpp:170
Definition: Tags.hpp:59
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Represents a multicolumn dat file inside an HDF5 file.
Definition: Dat.hpp:42
size_t get_l_max() const noexcept
retrieves the l_max that will be supplied to the DataBox in populate_hypersurface_boundary_data() ...
Definition: ReadBoundaryDataH5.hpp:589
bool populate_hypersurface_boundary_data(const gsl::not_null< db::DataBox< TagList > *> boundary_data_box, const double time) noexcept
Update the boundary_data_box entries for all tags in Tags::characteristic_worldtube_boundary_tags to ...
Definition: ReadBoundaryDataH5.hpp:416
Manages the &#39;reduced&#39; cached buffer dataset associated with a CCE worldtube and interpolates to reque...
Definition: ReadBoundaryDataH5.hpp:650
Create a span/view on a range, which is cheap to copy (one pointer).
Definition: Gsl.hpp:291
Prefix indicating a time derivative.
Definition: Prefixes.hpp:30
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
SpinWeighted< ComplexDataVector, Spin > inverse_swsh_transform(const size_t l_max, const size_t number_of_radial_points, const SpinWeighted< ComplexModalVector, Spin > &libsharp_coefficients) noexcept
Perform an inverse libsharp spin-weighted spherical harmonic transform on a single supplied SpinWeigh...
Definition: SwshTransform.cpp:118
bool time_is_outside_range(const double time) const noexcept override
The time can only be supported in the buffer update if it is between the first and last time of the i...
Definition: ReadBoundaryDataH5.hpp:307
Definition: ReadBoundaryDataH5.hpp:248
Definition: Determinant.hpp:11
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition...
Definition: SpinWeighted.hpp:25
constexpr size_t goldberg_mode_index(const size_t l_max, const size_t l, const int m, const size_t radial_offset=0) noexcept
Returns the index into a vector of modes consistent with .
Definition: SwshCoefficients.hpp:478
DataVector & get_time_buffer() noexcept override
The time buffer is supplied by non-const reference to allow views to easily point into the buffer...
Definition: ReadBoundaryDataH5.hpp:327
Definition: Tags.hpp:28
Abstract base class for utilities that are able to perform the buffer updating procedure needed by th...
Definition: ReadBoundaryDataH5.hpp:138
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Definition: Tags.hpp:54
Defines class h5::Version for storing version history of files.
Definition: DataBoxTag.hpp:29
double get_extraction_radius() const noexcept override
retrieves the extraction radius encoded in the filename
Definition: ReadBoundaryDataH5.hpp:204
void pup(PUP::er &p) noexcept
Serialization for Charm++.
Definition: ReadBoundaryDataH5.hpp:598
std::pair< size_t, size_t > get_time_span() const noexcept
retrieves the current time span associated with the buffer_updater_ for diagnostics ...
Definition: ReadBoundaryDataH5.hpp:687
Defines macro ASSERT.
double get_extraction_radius() const noexcept override
retrieves the extraction radius encoded in the filename
Definition: ReadBoundaryDataH5.hpp:316
#define WRAPPED_PUPable_abstract(className)
Wraps the Charm++ macro, see the Charm++ documentation.
Definition: CharmPupable.hpp:39
bool populate_hypersurface_boundary_data(gsl::not_null< db::DataBox< TagList > *> boundary_data_box, double time) noexcept
Update the boundary_data_box entries for all tags in Tags::characteristic_worldtube_boundary_tags to ...
std::pair< size_t, size_t > get_time_span() const noexcept
retrieves the current time span associated with the buffer_updater_ for diagnostics ...
Definition: ReadBoundaryDataH5.hpp:593
Stores a collection of function values.
Definition: DataVector.hpp:42
size_t get_l_max() const noexcept override
retrieves the l_max of the input file
Definition: ReadBoundaryDataH5.hpp:313
size_t get_l_max() const noexcept override
retrieves the l_max of the input file
Definition: ReadBoundaryDataH5.hpp:201
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:475
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
T read_data(const hid_t group_id, const std::string &dataset_name) noexcept
Read an array of rank 0-3 into an object.
Definition: Helpers.cpp:452
Defines class h5::Dat.
Marks an item as being a prefix to another tag.
Definition: DataBoxTag.hpp:111
void create_bondi_boundary_data_from_unnormalized_spec_modes(const gsl::not_null< db::DataBox< DataBoxTagList > *> bondi_boundary_data, const tnsr::ii< ComplexModalVector, 3 > &spatial_metric_coefficients, const tnsr::ii< ComplexModalVector, 3 > &dt_spatial_metric_coefficients, const tnsr::ii< ComplexModalVector, 3 > &dr_spatial_metric_coefficients, const tnsr::I< ComplexModalVector, 3 > &shift_coefficients, const tnsr::I< ComplexModalVector, 3 > &dt_shift_coefficients, const tnsr::I< ComplexModalVector, 3 > &dr_shift_coefficients, const Scalar< ComplexModalVector > &lapse_coefficients, const Scalar< ComplexModalVector > &dt_lapse_coefficients, const Scalar< ComplexModalVector > &dr_lapse_coefficients, const double extraction_radius, const size_t l_max) noexcept
Process the worldtube data from modal metric components and derivatives with incorrectly normalized r...
Definition: SpecBoundaryData.hpp:153
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Manages the cached buffer data associated with a CCE worldtube and interpolates to requested time poi...
Definition: ReadBoundaryDataH5.hpp:364
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182
A WorldtubeBufferUpdater specialized to the CCE input worldtube H5 file produced by SpEC...
Definition: ReadBoundaryDataH5.hpp:165
DataVector & get_time_buffer() noexcept override
The time buffer is supplied by non-const reference to allow views to easily point into the buffer...
Definition: ReadBoundaryDataH5.hpp:215