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