ReadSpecThirdOrderPiecewisePolynomial.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <memory>
9 #include <tuple>
10 #include <unordered_map>
11 #include <utility>
12 #include <vector>
13 
16 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
17 #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
18 #include "Domain/FunctionsOfTime/Tags.hpp"
19 #include "ErrorHandling/Assert.hpp"
20 #include "ErrorHandling/Error.hpp"
21 #include "IO/H5/AccessType.hpp"
22 #include "IO/H5/Dat.hpp"
23 #include "IO/H5/File.hpp"
24 #include "IO/Importers/Tags.hpp"
26 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
27 #include "Utilities/Gsl.hpp"
28 #include "Utilities/Requires.hpp"
29 #include "Utilities/TMPL.hpp"
30 #include "Utilities/TaggedTuple.hpp"
31 
32 namespace importers {
33 namespace Actions {
34 /// \brief Import SpEC `FunctionOfTime` data from an H5 file.
35 ///
36 /// Uses:
37 /// - DataBox:
38 /// - `importers::Tags::FunctionOfTimeFile`
39 /// - `importers::Tags::FunctionOfTimeNameMap`
40 ///
41 /// DataBox changes:
42 /// - Adds: nothing
43 /// - Removes: nothing
44 /// - Modifies:
45 /// - `::domain::Tags::FunctionsOfTime`
46 ///
47 /// Columns in the file to be read must have the following form:
48 /// - 0 = time
49 /// - 1 = time of last update
50 /// - 2 = number of components
51 /// - 3 = maximum derivative order
52 /// - 4 = version
53 /// - 5 = function
54 /// - 6 = d/dt (function)
55 /// - 7 = d^2/dt^2 (function)
56 /// - 8 = d^3/dt^3 (function)
57 ///
58 /// If the function has more than one component, columns 5-8 give
59 /// the first component and its derivatives, columns 9-12 give the second
60 /// component and its derivatives, etc.
61 ///
63  using const_global_cache_tags =
66  template <
67  typename DbTagsList, typename... InboxTags, typename Metavariables,
68  typename ArrayIndex, typename ActionList, typename ParallelComponent,
69  Requires<db::tag_is_retrievable_v<importers::Tags::FunctionOfTimeFile,
71  db::tag_is_retrievable_v<importers::Tags::FunctionOfTimeNameMap,
73  db::tag_is_retrievable_v<::domain::Tags::FunctionsOfTime,
74  db::DataBox<DbTagsList>>> = nullptr>
75  static auto apply(db::DataBox<DbTagsList>& box,
76  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
78  const ArrayIndex& /*array_index*/,
79  const ActionList /*meta*/,
80  const ParallelComponent* const /*meta*/) noexcept {
81  const std::string& file_name{
82  db::get<importers::Tags::FunctionOfTimeFile>(box)};
83 
84  // Get the map of SpEC -> SpECTRE FunctionofTime names
85  const std::map<std::string, std::string>& dataset_name_map{
86  db::get<importers::Tags::FunctionOfTimeNameMap>(box)};
87 
88  // Currently, only support order 3 piecewise polynomials.
89  // This could be generalized later, but the SpEC functions of time
90  // that we will read in with this action will always be 3rd-order
91  // piecewise polynomials
92  constexpr size_t max_deriv{3};
93 
96  spec_functions_of_time;
97 
99  std::make_unique<h5::H5File<h5::AccessType::ReadOnly>>(file_name);
100  for (const auto& spec_and_spectre_names : dataset_name_map) {
101  const std::string& spec_name = std::get<0>(spec_and_spectre_names);
102  const std::string& spectre_name = std::get<1>(spec_and_spectre_names);
103  // clang-tidy: use auto when initializing with a template cast to avoid
104  // duplicating the type name
105  const h5::Dat& dat_file = file->get<h5::Dat>("/" + spec_name); // NOLINT
106  const Matrix& dat_data = dat_file.get_data();
107 
108  // Check that the data in the file uses deriv order 3
109  // Column 3 of the file contains the derivative order
110  const size_t dat_max_deriv = dat_data(0, 3);
111  if (dat_max_deriv != max_deriv) {
112  file.reset();
113  ERROR("Deriv order in " << file_name << " should be " << max_deriv
114  << ", not " << dat_max_deriv);
115  }
116 
117  // Get the initial time ('time of last update') from the file
118  // and the values of the function and its derivatives at that time
119  const double start_time = dat_data(0, 1);
120 
121  // Currently, assume the same number of components are used
122  // at each time. This could be generalized if needed
123  const size_t number_of_components = dat_data(0, 2);
124 
125  std::array<DataVector, max_deriv + 1> initial_coefficients;
126  for (size_t deriv_order = 0; deriv_order < max_deriv + 1; ++deriv_order) {
127  gsl::at(initial_coefficients, deriv_order) =
128  DataVector(number_of_components);
129  for (size_t component = 0; component < number_of_components;
130  ++component) {
131  gsl::at(initial_coefficients, deriv_order)[component] =
132  dat_data(0, 5 + (max_deriv + 1) * component + deriv_order);
133  }
134  }
135  spec_functions_of_time[spectre_name] =
137  initial_coefficients);
138 
139  // Loop over the remaining times, updating the function of time
140  DataVector highest_derivative(number_of_components);
141  double time_last_updated = start_time;
142  for (size_t row = 1; row < dat_data.rows(); ++row) {
143  // If time of last update has changed, then update the FunctionOfTime
144  // The time of last update is stored in column 1 in the dat file
145  if (dat_data(row, 1) > time_last_updated) {
146  time_last_updated = dat_data(row, 1);
147  for (size_t a = 0; a < number_of_components; ++a) {
148  highest_derivative[a] =
149  dat_data(row, 5 + (max_deriv + 1) * a + max_deriv);
150  }
151  spec_functions_of_time[spectre_name].update(time_last_updated,
152  highest_derivative);
153  } else {
154  file.reset();
155  ERROR("Non-monotonic time found in FunctionOfTime data. "
156  << "Time " << dat_data(row, 1) << " follows time "
157  << time_last_updated << " while reading " << spectre_name
158  << "\n");
159  }
160  }
161  }
162 
163  // Mutate ::domain::Tags::FunctionsOfTime, adding the imported
164  // FunctionsOfTime to it
165  db::mutate<::domain::Tags::FunctionsOfTime>(
166  make_not_null(&box),
167  [&dataset_name_map, &spec_functions_of_time,
168  &file](const gsl::not_null<std::unordered_map<
169  std::string,
171  functions_of_time) {
172  for (auto spec_and_spectre_names : dataset_name_map) {
173  const std::string& spectre_name =
174  std::get<1>(spec_and_spectre_names);
175 
176  // The FunctionsOfTime we are mutating must already have
177  // an element with key==spectre_name; this action only
178  // mutates the value associated with that key
179  if ((*functions_of_time).count(spectre_name) == 0) {
180  file.reset();
181  ERROR("Trying to import data for key "
182  << spectre_name
183  << "in FunctionsOfTime, but FunctionsOfTime does not "
184  "contain that key.\n");
185  }
186  auto* piecewise_polynomial = dynamic_cast<
188  (*functions_of_time)[spectre_name].get());
189  if (piecewise_polynomial == nullptr) {
190  file.reset();
191  ERROR(
192  "The function of time with name spectre_name is not a "
193  "PiecewisePolynomial<"
194  << max_deriv
195  << "> and so cannot be set using "
196  "ReadSpecThirdOrderPiecewisePolynomial\n");
197  }
198  *piecewise_polynomial = spec_functions_of_time[spectre_name];
199  }
200  });
201  return std::forward_as_tuple(std::move(box));
202  }
203 };
204 } // namespace Actions
205 } // namespace importers
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
AccessType.hpp
std::string
utility
domain::Tags::FunctionsOfTime
The functions of time.
Definition: Tags.hpp:37
domain::FunctionsOfTime::PiecewisePolynomial
A function that has a piecewise-constant MaxDerivth derivative.
Definition: PiecewisePolynomial.hpp:22
h5::Dat
Represents a multicolumn dat file inside an HDF5 file.
Definition: Dat.hpp:42
vector
Error.hpp
importers::Tags::FunctionOfTimeFile
Path to an H5 file containing SpEC FunctionOfTime data to read.
Definition: Tags.hpp:176
importers
Items related to loading data from files.
Definition: ElementActions.hpp:18
tuple
File.hpp
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
DataBox.hpp
cstddef
Assert.hpp
array
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
std::map
memory
importers::Actions::ReadSpecThirdOrderPiecewisePolynomial
Import SpEC FunctionOfTime data from an H5 file.
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:62
importers::Tags::FunctionOfTimeNameMap
Pairs of strings mapping SpEC -> SpECTRE FunctionOfTime names.
Definition: Tags.hpp:196
Gsl.hpp
Dat.hpp
Requires.hpp
Matrix.hpp
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
h5::Dat::get_data
Matrix get_data() const
Definition: Dat.cpp:209
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
std::unique_ptr
unordered_map
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
TMPL.hpp
ConstGlobalCache.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183