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 : template <typename ResultDataType, size_t Dim, typename Frame>
74 1 : void interpolate_to_points(
75 : gsl::not_null<std::vector<ResultDataType>*> result,
76 : const std::variant<std::vector<std::string>, std::string>&
77 : volume_files_or_glob,
78 : const std::string& subfile_name, const ObservationVariant& observation,
79 : const std::vector<std::string>& tensor_components,
80 : const tnsr::I<DataVector, Dim, Frame>& target_points,
81 : bool extrapolate_into_excisions = false,
82 : bool error_on_missing_points = false,
83 : std::optional<size_t> num_threads = std::nullopt);
84 :
85 : template <size_t Dim, typename Frame>
86 1 : std::vector<DataVector> interpolate_to_points(
87 : const std::variant<std::vector<std::string>, std::string>&
88 : volume_files_or_glob,
89 : const std::string& subfile_name, const ObservationVariant& observation,
90 : const std::vector<std::string>& tensor_components,
91 : const tnsr::I<DataVector, Dim, Frame>& target_points,
92 : bool extrapolate_into_excisions = false,
93 : bool error_on_missing_points = false,
94 : std::optional<size_t> num_threads = std::nullopt) {
95 : std::vector<DataVector> interpolated_data{};
96 : interpolate_to_points(make_not_null(&interpolated_data), volume_files_or_glob,
97 : subfile_name, observation, tensor_components,
98 : target_points, extrapolate_into_excisions,
99 : error_on_missing_points, num_threads);
100 : return interpolated_data;
101 : }
102 :
103 : template <typename Tags, size_t Dim, typename Frame>
104 1 : tuples::tagged_tuple_from_typelist<Tags> interpolate_to_points(
105 : const std::variant<std::vector<std::string>, std::string>&
106 : volume_files_or_glob,
107 : const std::string& subfile_name, const ObservationVariant& observation,
108 : const tnsr::I<DataVector, Dim, Frame>& target_points,
109 : bool extrapolate_into_excisions = false,
110 : const bool error_on_missing_points = false,
111 : std::optional<size_t> num_threads = std::nullopt) {
112 : return make_tagged_tuple<Tags>(interpolate_to_points(
113 : volume_files_or_glob, subfile_name, observation,
114 : get_tensor_components<Tags>(), target_points, extrapolate_into_excisions,
115 : error_on_missing_points, num_threads));
116 : }
117 : /// @}
118 :
119 : /*!
120 : * \brief Interpolates data in volume files to target points by reading the
121 : * volume data into memory
122 : *
123 : * This class reads the volume data at the requested time into memory and can
124 : * then interpolate it to any number of target points.
125 : *
126 : * \par Thread safety
127 : * Constructing the `PointwiseInterpolator` on multiple threads at once is not
128 : * thread safe (unless HDF5 was built with thread-safety support) because the
129 : * constructor opens the H5 files and reads in the data. However, once the data
130 : * is loaded, the `interpolate_to_points` and `interpolate_to_point` functions
131 : * are thread safe. This means the volume data can be loaded once in a single
132 : * thread and then used by multiple threads to interpolate to different points
133 : * in parallel.
134 : */
135 : template <size_t Dim, typename Frame = ::Frame::Inertial>
136 1 : struct PointwiseInterpolator {
137 0 : PointwiseInterpolator() = default;
138 0 : PointwiseInterpolator(const std::variant<std::vector<std::string>,
139 : std::string>& volume_files_or_glob,
140 : const std::string& subfile_name,
141 : const ObservationVariant& observation,
142 : const std::vector<std::string>& tensor_components);
143 :
144 : /*!
145 : * \brief Interpolate to many points
146 : *
147 : * \param result the interpolated data at the target points. The outer vector
148 : * is the number of components and the inner vector is the number of target
149 : * points. Will be resized automatically.
150 : * \param target_points the points to interpolate to
151 : * \param extrapolate_into_excisions whether to extrapolate into excised
152 : * regions
153 : * \param error_on_missing_points whether to throw an error if any of the
154 : * target points are outside the domain
155 : * \param num_threads the number of OpenMP threads to use to parallelize the
156 : * interpolation over the target points. If not provided, the default number
157 : * of threads will be used.
158 : */
159 1 : void interpolate_to_points(
160 : gsl::not_null<std::vector<DataVector>*> result,
161 : const tnsr::I<DataVector, Dim, Frame>& target_points,
162 : bool extrapolate_into_excisions = false,
163 : bool error_on_missing_points = false,
164 : std::optional<size_t> num_threads = std::nullopt) const;
165 :
166 : /*!
167 : * \brief Interpolate to a single point
168 : *
169 : * This function is thread safe, so the interpolator can be constructed to
170 : * load the volume data once and then used by multiple threads to interpolate
171 : * to different points in parallel. Note that updating the `block_order` (if
172 : * provided) is not thread safe, so each thread should manage a separate
173 : * `block_order`.
174 : *
175 : * \param result the interpolated data at the target point. The vector is over
176 : * the number of components. Will be resized automatically.
177 : * \param target_point the point to interpolate to
178 : * \param block_order an optional priority order to search for the block
179 : * containing the target point. Will be updated when the point is found.
180 : * See `block_logical_coordinates_single_point` for more details.
181 : */
182 1 : void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
183 : const tnsr::I<double, Dim, Frame>& target_point,
184 : std::optional<gsl::not_null<std::vector<size_t>*>>
185 : block_order = std::nullopt) const;
186 :
187 : /*!
188 : * \brief Interpolate to a single point in block-logical coordinates
189 : *
190 : * \param result the interpolated data at the target point. The vector is over
191 : * the number of components. Will be resized automatically.
192 : * \param target_point the point to interpolate to in block-logical
193 : * coordinates.
194 : */
195 1 : void interpolate_to_point(
196 : gsl::not_null<std::vector<double>*> result,
197 : const IdPair<domain::BlockId,
198 : tnsr::I<double, Dim, ::Frame::BlockLogical>>& target_point)
199 : const;
200 :
201 0 : size_t obs_id() const { return obs_id_; }
202 0 : double time() const { return time_; }
203 0 : const Domain<Dim>& domain() const { return domain_; }
204 0 : const domain::FunctionsOfTimeMap& functions_of_time() const {
205 : return functions_of_time_;
206 : }
207 :
208 : private:
209 : // Loaded from data files
210 0 : size_t obs_id_ = std::numeric_limits<size_t>::max();
211 0 : double time_ = std::numeric_limits<double>::signaling_NaN();
212 0 : Domain<Dim> domain_;
213 0 : domain::FunctionsOfTimeMap functions_of_time_;
214 : // Outer vector is the source data file
215 0 : std::vector<std::vector<ElementId<Dim>>> element_ids_;
216 : std::vector<std::map<size_t, domain::ElementSearchTree<Dim>>>
217 0 : element_search_trees_;
218 : std::vector<
219 : std::unordered_map<ElementId<Dim>, std::tuple<Mesh<Dim>, size_t, size_t>>>
220 0 : meshes_;
221 : std::vector<std::vector<std::variant<DataVector, std::vector<float>>>>
222 0 : tensor_data_;
223 : };
224 :
225 : } // namespace spectre::Exporter
|
