Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <string>
10 : #include <unordered_map>
11 : #include <variant>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/TagName.hpp"
15 : #include "DataStructures/TaggedTuple.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "Domain/BlockLogicalCoordinates.hpp"
18 : #include "Domain/Domain.hpp"
19 : #include "Domain/Structure/ElementSearchTree.hpp"
20 : #include "IO/Exporter/Exporter.hpp"
21 :
22 : namespace spectre::Exporter {
23 :
24 : /// Collect tensor component names from Tags list
25 : template <typename Tags>
26 1 : auto get_tensor_components() {
27 : std::vector<std::string> tensor_components{};
28 : tmpl::for_each<Tags>([&tensor_components](auto tag_v) {
29 : using tensor_tag = tmpl::type_from<decltype(tag_v)>;
30 : using TensorType = typename tensor_tag::type;
31 : const std::string tag_name = db::tag_name<tensor_tag>();
32 : for (size_t i = 0; i < TensorType::size(); ++i) {
33 : const std::string component_name =
34 : tag_name + TensorType::component_suffix(i);
35 : tensor_components.push_back(component_name);
36 : }
37 : });
38 : return tensor_components;
39 : }
40 :
41 : /// Convert tensor components to a tagged_tuple
42 : template <typename Tags, typename DataType>
43 1 : auto make_tagged_tuple(std::vector<DataType> interpolated_data) {
44 : tuples::tagged_tuple_from_typelist<Tags> result{};
45 : size_t component_index = 0;
46 : tmpl::for_each<Tags>(
47 : [&component_index, &interpolated_data, &result](auto tag_v) {
48 : using tensor_tag = tmpl::type_from<decltype(tag_v)>;
49 : using TensorType = typename tensor_tag::type;
50 : for (size_t i = 0; i < TensorType::size(); ++i) {
51 : get<tensor_tag>(result)[i] =
52 : std::move(interpolated_data[component_index]);
53 : ++component_index;
54 : }
55 : });
56 : return result;
57 : }
58 :
59 : /// @{
60 : /*!
61 : * \brief Interpolates data in volume files to target points
62 : *
63 : * These are overloads of the `interpolate_to_points` function that work with
64 : * Tensor types and tags, rather than the raw C++ types that are used in
65 : * Exporter.hpp so it can be used by external programs.
66 : *
67 : * The `Tags` template parameter is a typelist of tags that should be read from
68 : * the volume files. The dataset names to read are constructed from the tag
69 : * names. Here is an example of how to use this function:
70 : *
71 : * \snippet Test_Exporter.cpp interpolate_tensors_to_points_example
72 : *
73 : * The overload that takes the `result` as an argument resizes it to the number
74 : * of tensor components and target points. If it already has the correct size,
75 : * the data is written into it directly, so its components can be non-owning
76 : * `DataVector`s that reference external memory.
77 : */
78 : template <typename ResultDataType, size_t Dim, typename Frame>
79 1 : void interpolate_to_points(
80 : gsl::not_null<std::vector<ResultDataType>*> result,
81 : const std::variant<std::vector<std::string>, std::string>&
82 : volume_files_or_glob,
83 : const std::string& subfile_name, const ObservationVariant& observation,
84 : const std::vector<std::string>& tensor_components,
85 : const tnsr::I<DataVector, Dim, Frame>& target_points,
86 : bool extrapolate_into_excisions = false,
87 : bool error_on_missing_points = false,
88 : std::optional<size_t> num_threads = std::nullopt);
89 :
90 : template <size_t Dim, typename Frame>
91 1 : std::vector<DataVector> interpolate_to_points(
92 : const std::variant<std::vector<std::string>, std::string>&
93 : volume_files_or_glob,
94 : const std::string& subfile_name, const ObservationVariant& observation,
95 : const std::vector<std::string>& tensor_components,
96 : const tnsr::I<DataVector, Dim, Frame>& target_points,
97 : bool extrapolate_into_excisions = false,
98 : bool error_on_missing_points = false,
99 : std::optional<size_t> num_threads = std::nullopt) {
100 : std::vector<DataVector> interpolated_data{};
101 : interpolate_to_points(make_not_null(&interpolated_data), volume_files_or_glob,
102 : subfile_name, observation, tensor_components,
103 : target_points, extrapolate_into_excisions,
104 : error_on_missing_points, num_threads);
105 : return interpolated_data;
106 : }
107 :
108 : template <typename Tags, size_t Dim, typename Frame>
109 1 : tuples::tagged_tuple_from_typelist<Tags> interpolate_to_points(
110 : const std::variant<std::vector<std::string>, std::string>&
111 : volume_files_or_glob,
112 : const std::string& subfile_name, const ObservationVariant& observation,
113 : const tnsr::I<DataVector, Dim, Frame>& target_points,
114 : bool extrapolate_into_excisions = false,
115 : const bool error_on_missing_points = false,
116 : std::optional<size_t> num_threads = std::nullopt) {
117 : return make_tagged_tuple<Tags>(interpolate_to_points(
118 : volume_files_or_glob, subfile_name, observation,
119 : get_tensor_components<Tags>(), target_points, extrapolate_into_excisions,
120 : error_on_missing_points, num_threads));
121 : }
122 : /// @}
123 :
124 : /*!
125 : * \brief Interpolates data in volume files to target points by reading the
126 : * volume data into memory
127 : *
128 : * This class reads the volume data at the requested time into memory and can
129 : * then interpolate it to any number of target points.
130 : *
131 : * \par Thread safety
132 : * Constructing the `PointwiseInterpolator` on multiple threads at once is not
133 : * thread safe (unless HDF5 was built with thread-safety support) because the
134 : * constructor opens the H5 files and reads in the data. However, once the data
135 : * is loaded, the `interpolate_to_points` and `interpolate_to_point` functions
136 : * are thread safe. This means the volume data can be loaded once in a single
137 : * thread and then used by multiple threads to interpolate to different points
138 : * in parallel.
139 : */
140 : template <size_t Dim, typename Frame = ::Frame::Inertial>
141 1 : struct PointwiseInterpolator {
142 0 : PointwiseInterpolator() = default;
143 0 : PointwiseInterpolator(const std::variant<std::vector<std::string>,
144 : std::string>& volume_files_or_glob,
145 : const std::string& subfile_name,
146 : const ObservationVariant& observation,
147 : const std::vector<std::string>& tensor_components);
148 :
149 : /*!
150 : * \brief Interpolate to many points
151 : *
152 : * \param result the interpolated data at the target points. The outer vector
153 : * is the number of components and the inner vector is the number of target
154 : * points. Will be resized automatically.
155 : * \param target_points the points to interpolate to
156 : * \param extrapolate_into_excisions whether to extrapolate into excised
157 : * regions
158 : * \param error_on_missing_points whether to throw an error if any of the
159 : * target points are outside the domain
160 : * \param num_threads the number of OpenMP threads to use to parallelize the
161 : * interpolation over the target points. If not provided, the default number
162 : * of threads will be used.
163 : */
164 1 : void interpolate_to_points(
165 : gsl::not_null<std::vector<DataVector>*> result,
166 : const tnsr::I<DataVector, Dim, Frame>& target_points,
167 : bool extrapolate_into_excisions = false,
168 : bool error_on_missing_points = false,
169 : std::optional<size_t> num_threads = std::nullopt) const;
170 :
171 : /*!
172 : * \brief Interpolate to a single point
173 : *
174 : * This function is thread safe, so the interpolator can be constructed to
175 : * load the volume data once and then used by multiple threads to interpolate
176 : * to different points in parallel. Note that updating the `block_order` (if
177 : * provided) is not thread safe, so each thread should manage a separate
178 : * `block_order`.
179 : *
180 : * \param result the interpolated data at the target point. The vector is over
181 : * the number of components. Will be resized automatically.
182 : * \param target_point the point to interpolate to
183 : * \param block_order an optional priority order to search for the block
184 : * containing the target point. Will be updated when the point is found.
185 : * See `block_logical_coordinates_single_point` for more details.
186 : */
187 1 : void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
188 : const tnsr::I<double, Dim, Frame>& target_point,
189 : std::optional<gsl::not_null<std::vector<size_t>*>>
190 : block_order = std::nullopt) const;
191 :
192 : /*!
193 : * \brief Interpolate to a single point in block-logical coordinates
194 : *
195 : * \param result the interpolated data at the target point. The vector is over
196 : * the number of components. Will be resized automatically.
197 : * \param target_point the point to interpolate to in block-logical
198 : * coordinates.
199 : */
200 1 : void interpolate_to_point(
201 : gsl::not_null<std::vector<double>*> result,
202 : const IdPair<domain::BlockId,
203 : tnsr::I<double, Dim, ::Frame::BlockLogical>>& target_point)
204 : const;
205 :
206 0 : size_t obs_id() const { return obs_id_; }
207 0 : double time() const { return time_; }
208 0 : const Domain<Dim>& domain() const { return domain_; }
209 0 : const domain::FunctionsOfTimeMap& functions_of_time() const {
210 : return functions_of_time_;
211 : }
212 :
213 : private:
214 : // Loaded from data files
215 0 : size_t obs_id_ = std::numeric_limits<size_t>::max();
216 0 : double time_ = std::numeric_limits<double>::signaling_NaN();
217 0 : Domain<Dim> domain_;
218 0 : domain::FunctionsOfTimeMap functions_of_time_;
219 : // Outer vector is the source data file
220 0 : std::vector<std::vector<ElementId<Dim>>> element_ids_;
221 : std::vector<std::map<size_t, domain::ElementSearchTree<Dim>>>
222 0 : element_search_trees_;
223 : std::vector<
224 : std::unordered_map<ElementId<Dim>, std::tuple<Mesh<Dim>, size_t, size_t>>>
225 0 : meshes_;
226 : std::vector<std::vector<std::variant<DataVector, std::vector<float>>>>
227 0 : tensor_data_;
228 : };
229 :
230 : } // namespace spectre::Exporter
|
