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