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 <limits>
9 : #include <optional>
10 : #include <pup.h>
11 : #include <set>
12 : #include <sstream>
13 : #include <string>
14 : #include <utility>
15 : #include <vector>
16 :
17 : #include "DataStructures/DataBox/DataBox.hpp"
18 : #include "DataStructures/DataVector.hpp"
19 : #include "DataStructures/Index.hpp"
20 : #include "Domain/BlockLogicalCoordinates.hpp"
21 : #include "Domain/Creators/Tags/Domain.hpp"
22 : #include "Domain/Domain.hpp"
23 : #include "Domain/ElementLogicalCoordinates.hpp"
24 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
25 : #include "Domain/Structure/BlockId.hpp"
26 : #include "Domain/Structure/ElementId.hpp"
27 : #include "Domain/Tags.hpp"
28 : #include "IO/Logging/Verbosity.hpp"
29 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
30 : #include "Options/String.hpp"
31 : #include "Parallel/GlobalCache.hpp"
32 : #include "Parallel/Invoke.hpp"
33 : #include "Parallel/ParallelComponentHelpers.hpp"
34 : #include "Parallel/Printf/Printf.hpp"
35 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
36 : #include "ParallelAlgorithms/Interpolation/Actions/InterpolationTargetVarsFromElement.hpp"
37 : #include "ParallelAlgorithms/Interpolation/Events/GetComputeItemsOnSource.hpp"
38 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
39 : #include "ParallelAlgorithms/Interpolation/PointInfoTag.hpp"
40 : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
41 : #include "ParallelAlgorithms/Interpolation/Targets/Sphere.hpp"
42 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
43 : #include "Utilities/Algorithm.hpp"
44 : #include "Utilities/CartesianProduct.hpp"
45 : #include "Utilities/ConstantExpressions.hpp"
46 : #include "Utilities/PrettyType.hpp"
47 : #include "Utilities/Serialization/CharmPupable.hpp"
48 : #include "Utilities/TMPL.hpp"
49 : #include "Utilities/TypeTraits/CreateGetStaticMemberVariableOrDefault.hpp"
50 : #include "Utilities/TypeTraits/IsA.hpp"
51 :
52 : /// \cond
53 : namespace Events::Tags {
54 : template <size_t Dim>
55 : struct ObserverMesh;
56 : template <size_t Dim, typename Fr>
57 : struct ObserverCoordinates;
58 : } // namespace Events::Tags
59 : namespace Tags {
60 : template <typename TagsList>
61 : struct Variables;
62 : } // namespace Tags
63 : template <size_t VolumeDim>
64 : class ElementId;
65 : namespace intrp {
66 : template <typename Metavariables, typename Tag>
67 : struct InterpolationTarget;
68 : } // namespace intrp
69 : /// \endcond
70 :
71 : namespace intrp::Events {
72 : /// \cond
73 : template <size_t VolumeDim, typename InterpolationTargetTag,
74 : typename SourceVarTags>
75 : class InterpolateWithoutInterpComponent;
76 : /// \endcond
77 :
78 : /*!
79 : * \brief Does an interpolation onto an InterpolationTargetTag by calling
80 : * Actions on the InterpolationTarget component.
81 : *
82 : * \note The `intrp::TargetPoints::Sphere` target is handled specially because
83 : * it has the potential to be very slow due to it usually having the most points
84 : * out of all the stationary targets. An optimization for the future would be to
85 : * have each target be responsible for intelligently computing the
86 : * `block_logical_coordinates` for it's own points.
87 : */
88 : template <size_t VolumeDim, typename InterpolationTargetTag,
89 : typename... SourceVarTags>
90 1 : class InterpolateWithoutInterpComponent<VolumeDim, InterpolationTargetTag,
91 : tmpl::list<SourceVarTags...>>
92 : : public Event {
93 : private:
94 0 : using frame = typename InterpolationTargetTag::compute_target_points::frame;
95 :
96 : public:
97 : /// \cond
98 : explicit InterpolateWithoutInterpComponent(CkMigrateMessage* /*unused*/) {}
99 : using PUP::able::register_constructor;
100 : WRAPPED_PUPable_decl_template(InterpolateWithoutInterpComponent); // NOLINT
101 : /// \endcond
102 :
103 0 : using options = tmpl::list<>;
104 0 : static constexpr Options::String help =
105 : "Does interpolation using the given InterpolationTargetTag, "
106 : "without an Interpolator ParallelComponent.";
107 :
108 0 : static std::string name() {
109 : return pretty_type::name<InterpolationTargetTag>();
110 : }
111 :
112 0 : InterpolateWithoutInterpComponent() = default;
113 :
114 0 : using compute_tags_for_observation_box =
115 : detail::get_compute_items_on_source_or_default_t<InterpolationTargetTag,
116 : tmpl::list<>>;
117 :
118 0 : using return_tags = tmpl::list<>;
119 0 : using argument_tags = tmpl::list<
120 : typename InterpolationTargetTag::temporal_id,
121 : Tags::PointInfo<InterpolationTargetTag, tmpl::size_t<VolumeDim>>,
122 : ::Events::Tags::ObserverMesh<VolumeDim>,
123 : // We always grab the DG coords because we use them to create a
124 : // bounding box for the sphere target optimization. DG coords
125 : // have points on the boundary, while FD coords don't. If we
126 : // had used FD coords, it's possible a target point would fall
127 : // between the outermost gridpoints of two elements. This point
128 : // would be outside our bounding box, and thus wouldn't get
129 : // interpolated to. We avoid this by always using DG coords,
130 : // even if the mesh is FD.
131 : domain::Tags::Coordinates<VolumeDim, frame>, SourceVarTags...>;
132 :
133 : template <typename ParallelComponent, typename Metavariables>
134 0 : void operator()(
135 : const typename InterpolationTargetTag::temporal_id::type& temporal_id,
136 : const tnsr::I<DataVector, VolumeDim, frame>& all_target_points,
137 : const Mesh<VolumeDim>& mesh,
138 : const tnsr::I<DataVector, VolumeDim, frame> coordinates,
139 : const typename SourceVarTags::type&... source_vars_input,
140 : Parallel::GlobalCache<Metavariables>& cache,
141 : const ElementId<VolumeDim>& array_index,
142 : const ParallelComponent* const /*meta*/,
143 : const ObservationValue& /*observation_value*/) const {
144 : std::vector<BlockLogicalCoords<VolumeDim>> block_logical_coords{};
145 :
146 : std::stringstream ss{};
147 : ss << std::setprecision(std::numeric_limits<double>::digits10 + 4)
148 : << std::scientific;
149 : const ::Verbosity verbosity = Parallel::get<intrp::Tags::Verbosity>(cache);
150 : const bool debug_print = verbosity >= ::Verbosity::Debug;
151 : if (debug_print) {
152 : ss << InterpolationTarget_detail::target_output_prefix<
153 : InterpolateWithoutInterpComponent, InterpolationTargetTag>(
154 : temporal_id)
155 : << ", " << array_index << ", ";
156 : }
157 :
158 : // The sphere target is special because we have a better idea of where the
159 : // points will be.
160 : if constexpr (tt::is_a_v<
161 : TargetPoints::Sphere,
162 : typename InterpolationTargetTag::compute_target_points>) {
163 : static_assert(VolumeDim == 3,
164 : "Sphere target can only be used for VolumeDim = 3.");
165 : // Extremum for x, y, z, r
166 : std::array<std::pair<double, double>, 4> min_max_coordinates{};
167 :
168 : const auto& sphere =
169 : Parallel::get<Tags::Sphere<InterpolationTargetTag>>(cache);
170 : const std::array<double, 3>& center = sphere.center;
171 :
172 : // Calculate r^2 from center of sphere because sqrt is expensive
173 : DataVector radii_squared{get<0>(coordinates).size(), 0.0};
174 : for (size_t i = 0; i < VolumeDim; i++) {
175 : radii_squared += square(coordinates.get(i) - gsl::at(center, i));
176 : }
177 :
178 : // Compute min and max
179 : {
180 : const auto [min, max] = alg::minmax_element(radii_squared);
181 : min_max_coordinates[3].first = *min;
182 : min_max_coordinates[3].second = *max;
183 : }
184 :
185 : const std::set<double>& radii_of_sphere_target = sphere.radii;
186 : const size_t l_max = sphere.l_max;
187 :
188 : const size_t number_of_angular_points = (l_max + 1) * (2 * l_max + 1);
189 : // first size_t = position of first radius in bounds
190 : // second size_t = total bounds to use/check
191 : std::optional<std::pair<size_t, size_t>> offset_and_num_points{};
192 :
193 : // Have a very small buffer just in case of roundoff
194 : double epsilon =
195 : (min_max_coordinates[3].second - min_max_coordinates[3].first) *
196 : std::numeric_limits<double>::epsilon() * 100.0;
197 : size_t offset_index = 0;
198 : // Check if any radii of the target are within the radii of our element
199 : for (double radius : radii_of_sphere_target) {
200 : const double square_radius = square(radius);
201 : if (square_radius >=
202 : (gsl::at(min_max_coordinates, 3).first - epsilon) and
203 : square_radius <=
204 : (gsl::at(min_max_coordinates, 3).second + epsilon)) {
205 : if (offset_and_num_points.has_value()) {
206 : offset_and_num_points->second += number_of_angular_points;
207 : } else {
208 : offset_and_num_points =
209 : std::make_pair(offset_index * number_of_angular_points,
210 : number_of_angular_points);
211 : }
212 : }
213 : offset_index++;
214 : }
215 :
216 : // If no radii pass through this element, there's nothing to do so return
217 : if (not offset_and_num_points.has_value()) {
218 : if (debug_print) {
219 : using ::operator<<;
220 : ss << "No radii in this element.";
221 : Parallel::printf("%s\n", ss.str());
222 : }
223 : return;
224 : }
225 :
226 : // Get the x,y,z bounds
227 : for (size_t i = 0; i < VolumeDim; i++) {
228 : const auto [min, max] = alg::minmax_element(coordinates.get(i));
229 : gsl::at(min_max_coordinates, i).first = *min;
230 : gsl::at(min_max_coordinates, i).second = *max;
231 : }
232 :
233 : const tnsr::I<DataVector, VolumeDim, frame> target_points_to_check{};
234 : // Use the offset and number of points to create a view. We assume that if
235 : // there are multiple radii in this element, that they are successive
236 : // radii. If this wasn't true, that'd be a really weird topology.
237 : for (size_t i = 0; i < VolumeDim; i++) {
238 : make_const_view(make_not_null(&target_points_to_check.get(i)),
239 : all_target_points.get(i), offset_and_num_points->first,
240 : offset_and_num_points->second);
241 : }
242 :
243 : // To break out of inner loop and skip the point
244 : bool skip_point = false;
245 : tnsr::I<double, VolumeDim, frame> sphere_coords_to_map{};
246 : // Other code below expects block_logical_coords to be sized to the total
247 : // number of points on the target (including all radii). By default these
248 : // will all be nullopt and we'll only fill the ones we need
249 : block_logical_coords.resize(get<0>(all_target_points).size());
250 :
251 : const Block<VolumeDim>& block =
252 : Parallel::get<domain::Tags::Domain<VolumeDim>>(cache)
253 : .blocks()[array_index.block_id()];
254 :
255 : // Now for every radius in this element, we check if their points are
256 : // within the x,y,z bounds of the element. If they are, map the point to
257 : // the block logical frame and add it to the vector f all block logical
258 : // coordinates.
259 : for (size_t index = 0; index < get<0>(target_points_to_check).size();
260 : index++) {
261 : skip_point = false;
262 : for (size_t i = 0; i < VolumeDim; i++) {
263 : const double coord = target_points_to_check.get(i)[index];
264 : epsilon = (gsl::at(min_max_coordinates, i).second -
265 : gsl::at(min_max_coordinates, i).first) *
266 : std::numeric_limits<double>::epsilon() * 100.0;
267 : // If a point is outside any of the bounding box, skip it
268 : if (coord < (gsl::at(min_max_coordinates, i).first - epsilon) or
269 : coord > (gsl::at(min_max_coordinates, i).second + epsilon)) {
270 : skip_point = true;
271 : break;
272 : }
273 :
274 : sphere_coords_to_map.get(i) = coord;
275 : }
276 :
277 : if (skip_point) {
278 : continue;
279 : }
280 :
281 : std::optional<tnsr::I<double, VolumeDim, ::Frame::BlockLogical>>
282 : block_coords_of_target_point{};
283 :
284 : if constexpr (Parallel::is_in_global_cache<
285 : Metavariables, domain::Tags::FunctionsOfTime>) {
286 : const auto& functions_of_time =
287 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
288 : const double time =
289 : InterpolationTarget_detail::get_temporal_id_value(temporal_id);
290 :
291 : block_coords_of_target_point = block_logical_coordinates_single_point(
292 : sphere_coords_to_map, block, time, functions_of_time);
293 : } else {
294 : block_coords_of_target_point = block_logical_coordinates_single_point(
295 : sphere_coords_to_map, block);
296 : }
297 :
298 : if (block_coords_of_target_point.has_value()) {
299 : // Get index into vector of all grid points of the target. This is
300 : // just the offset + index
301 : block_logical_coords[offset_and_num_points->first + index] =
302 : make_id_pair(domain::BlockId(array_index.block_id()),
303 : std::move(block_coords_of_target_point.value()));
304 : }
305 : }
306 : } else {
307 : (void)coordinates;
308 : block_logical_coords = InterpolationTarget_detail::block_logical_coords<
309 : InterpolationTargetTag>(cache, all_target_points, temporal_id);
310 : }
311 :
312 : const std::vector<ElementId<VolumeDim>> element_ids{{array_index}};
313 : const auto element_coord_holders =
314 : element_logical_coordinates(element_ids, block_logical_coords);
315 :
316 : if (element_coord_holders.count(array_index) == 0) {
317 : if (debug_print) {
318 : ss << "No target points in this element. Skipping.";
319 : Parallel::printf("%s\n", ss.str());
320 : }
321 : // There are no target points in this element, so we don't need
322 : // to do anything.
323 : return;
324 : }
325 :
326 : // There are points in this element, so interpolate to them and
327 : // send the interpolated data to the target. This is done
328 : // in several steps:
329 : const auto& element_coord_holder = element_coord_holders.at(array_index);
330 :
331 : // 1. Get the list of variables
332 : Variables<typename InterpolationTargetTag::vars_to_interpolate_to_target>
333 : interp_vars(mesh.number_of_grid_points());
334 :
335 : if constexpr (InterpolationTarget_detail::has_compute_vars_to_interpolate_v<
336 : InterpolationTargetTag>) {
337 : // 1a. Call compute_vars_to_interpolate. Need the source in a
338 : // Variables, so copy the variables here.
339 : // This copy would be unnecessary if we passed a Variables into
340 : // InterpolateWithoutInterpComponent instead of passing
341 : // individual Tensors, which would require that this Variables is
342 : // something in the DataBox. (Note that
343 : // InterpolationTarget_detail::compute_dest_vars_from_source_vars
344 : // allows the source variables to be different from the
345 : // destination variables).
346 : Variables<tmpl::list<SourceVarTags...>> source_vars(
347 : mesh.number_of_grid_points());
348 : [[maybe_unused]] const auto copy_to_variables =
349 : [&source_vars](const auto source_var_tag_v, const auto& source_var) {
350 : using source_var_tag = tmpl::type_from<decltype(source_var_tag_v)>;
351 : get<source_var_tag>(source_vars) = source_var;
352 : return 0;
353 : };
354 : expand_pack(copy_to_variables(tmpl::type_<SourceVarTags>{},
355 : source_vars_input)...);
356 :
357 : InterpolationTarget_detail::compute_dest_vars_from_source_vars<
358 : InterpolationTargetTag>(make_not_null(&interp_vars), source_vars,
359 : get<domain::Tags::Domain<VolumeDim>>(cache),
360 : mesh, array_index, cache, temporal_id);
361 : } else {
362 : // 1b. There is no compute_vars_to_interpolate. So copy the
363 : // source vars directly into the variables.
364 : // This copy would be unnecessary if:
365 : // - We passed a Variables into InterpolateWithoutInterpComponent
366 : // instead of passing individual Tensors.
367 : // and
368 : // - This Variables was actually something in the DataBox.
369 : // and
370 : // - Either the passed-in Variables was exactly the same as
371 : // InterpolationTargetTag::vars_to_interpolate_to_target,
372 : // or IrregularInterpolant::interpolate had the ability to
373 : // interpolate only a subset of the Variables passed into it,
374 : // or IrregularInterpolant::interpolate can interpolate individual
375 : // DataVectors.
376 : [[maybe_unused]] const auto copy_to_variables =
377 : [&interp_vars](const auto tensor_tag_v, const auto& tensor) {
378 : using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
379 : get<tensor_tag>(interp_vars) = tensor;
380 : return 0;
381 : };
382 : expand_pack(copy_to_variables(tmpl::type_<SourceVarTags>{},
383 : source_vars_input)...);
384 : }
385 :
386 : // 2. Set up interpolator
387 : intrp::Irregular<VolumeDim> interpolator(
388 : mesh, element_coord_holder.element_logical_coords);
389 :
390 : // 3. Interpolate and send interpolated data to target
391 : auto& receiver_proxy = Parallel::get_parallel_component<
392 : InterpolationTarget<Metavariables, InterpolationTargetTag>>(cache);
393 : Parallel::simple_action<
394 : Actions::InterpolationTargetVarsFromElement<InterpolationTargetTag>>(
395 : receiver_proxy,
396 : std::vector<Variables<
397 : typename InterpolationTargetTag::vars_to_interpolate_to_target>>(
398 : {interpolator.interpolate(interp_vars)}),
399 : block_logical_coords,
400 : std::vector<std::vector<size_t>>({element_coord_holder.offsets}),
401 : temporal_id);
402 :
403 : if (debug_print) {
404 : ss << "Sending points and vars to target.";
405 : Parallel::printf("%s\n", ss.str());
406 : }
407 : }
408 :
409 0 : using is_ready_argument_tags = tmpl::list<>;
410 :
411 : template <typename ArrayIndex, typename Component, typename Metavariables>
412 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
413 : const ArrayIndex& /*array_index*/,
414 : const Component* const /*meta*/) const {
415 : return true;
416 : }
417 :
418 1 : bool needs_evolved_variables() const override { return true; }
419 : };
420 :
421 : /// \cond
422 : template <size_t VolumeDim, typename InterpolationTargetTag,
423 : typename... SourceVarTags>
424 : PUP::able::PUP_ID
425 : InterpolateWithoutInterpComponent<VolumeDim, InterpolationTargetTag,
426 : tmpl::list<SourceVarTags...>>::my_PUP_ID =
427 : 0; // NOLINT
428 : /// \endcond
429 :
430 : } // namespace intrp::Events
|
