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