Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <type_traits>
9 : #include <unordered_set>
10 : #include <utility>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/IdPair.hpp"
15 : #include "DataStructures/LinkedMessageId.hpp"
16 : #include "DataStructures/Tensor/Metafunctions.hpp"
17 : #include "DataStructures/VariablesTag.hpp"
18 : #include "Domain/BlockLogicalCoordinates.hpp"
19 : #include "Domain/CoordinateMaps/Composition.hpp"
20 : #include "Domain/Creators/Tags/Domain.hpp"
21 : #include "Domain/ElementToBlockLogicalMap.hpp"
22 : #include "Domain/TagsTimeDependent.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Parallel/Invoke.hpp"
25 : #include "Parallel/ParallelComponentHelpers.hpp"
26 : #include "ParallelAlgorithms/Interpolation/TagsMetafunctions.hpp"
27 : #include "Utilities/Algorithm.hpp"
28 : #include "Utilities/Gsl.hpp"
29 : #include "Utilities/Literals.hpp"
30 : #include "Utilities/Requires.hpp"
31 : #include "Utilities/TMPL.hpp"
32 : #include "Utilities/TaggedTuple.hpp"
33 : #include "Utilities/TypeTraits.hpp"
34 : #include "Utilities/TypeTraits/CreateHasStaticMemberVariable.hpp"
35 : #include "Utilities/TypeTraits/CreateHasTypeAlias.hpp"
36 : #include "Utilities/TypeTraits/CreateIsCallable.hpp"
37 : #include "Utilities/TypeTraits/IsA.hpp"
38 :
39 : /// \cond
40 : // IWYU pragma: no_forward_declare db::DataBox
41 : namespace intrp {
42 : template <class Metavariables>
43 : struct Interpolator;
44 : template <typename Metavariables, typename InterpolationTargetTag>
45 : class InterpolationTarget;
46 : namespace Actions {
47 : template <typename InterpolationTargetTag>
48 : struct CleanUpInterpolator;
49 : template <typename InterpolationTargetTag>
50 : struct ReceivePoints;
51 : } // namespace Actions
52 : namespace Tags {
53 : template <typename TemporalId>
54 : struct IndicesOfFilledInterpPoints;
55 : template <typename TemporalId>
56 : struct IndicesOfInvalidInterpPoints;
57 : template <typename InterpolationTargetTag, typename TemporalId>
58 : struct InterpolatedVars;
59 : template <typename TemporalId>
60 : struct CompletedTemporalIds;
61 : template <typename TemporalId>
62 : struct PendingTemporalIds;
63 : template <typename TemporalId>
64 : struct TemporalIds;
65 : } // namespace Tags
66 : namespace TargetPoints {
67 : template <typename InterpolationTargetTag, typename Frame>
68 : struct Sphere;
69 : } // namespace TargetPoints
70 : } // namespace intrp
71 : template <typename Id>
72 : struct LinkedMessageId;
73 : class TimeStepId;
74 : template <typename TagsList>
75 : struct Variables;
76 : /// \endcond
77 :
78 : namespace intrp {
79 :
80 : namespace InterpolationTarget_detail {
81 : CREATE_IS_CALLABLE(substep_time)
82 : CREATE_IS_CALLABLE_V(substep_time)
83 : template <typename T>
84 : double get_temporal_id_value(const T& id) {
85 : if constexpr (std::is_same_v<T, LinkedMessageId<double>>) {
86 : return id.id;
87 : } else if constexpr (is_substep_time_callable_v<T>) {
88 : return id.substep_time();
89 : } else {
90 : static_assert(std::is_same_v<T, double>,
91 : "get_temporal_id_value only supports 'double', "
92 : "'LinkedMessageId<double>', or 'TimeStepId'.");
93 : return id;
94 : }
95 : }
96 :
97 : CREATE_IS_CALLABLE(apply)
98 :
99 : template <typename T, typename DbTags, typename Metavariables,
100 : typename TemporalId>
101 : bool apply_callbacks(
102 : const gsl::not_null<db::DataBox<DbTags>*> box,
103 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
104 : const TemporalId& temporal_id) {
105 : using callbacks = typename T::post_interpolation_callbacks;
106 : constexpr bool all_simple_callbacks =
107 : tmpl::all<callbacks,
108 : is_apply_callable<tmpl::_1, tmpl::pin<decltype(*box)>,
109 : tmpl::pin<decltype(*cache)>,
110 : tmpl::pin<decltype(temporal_id)>>>::value;
111 : if constexpr (all_simple_callbacks) {
112 : tmpl::for_each<callbacks>([&](const auto callback_v) {
113 : using callback = tmpl::type_from<decltype(callback_v)>;
114 : callback::apply(*box, *cache, temporal_id);
115 : });
116 :
117 : // For the simpler callback function, we will always clean up volume data,
118 : // so we return true here.
119 : return true;
120 : } else {
121 : static_assert(tmpl::size<callbacks>::value == 1,
122 : "Can only have 1 post_interpolation_callbacks that mutates "
123 : "the target box.");
124 : return tmpl::front<callbacks>::apply(box, cache, temporal_id);
125 : }
126 : }
127 :
128 : CREATE_HAS_STATIC_MEMBER_VARIABLE(fill_invalid_points_with)
129 : CREATE_HAS_STATIC_MEMBER_VARIABLE_V(fill_invalid_points_with)
130 :
131 : template <typename T, bool B>
132 : struct fill_value_or_no_such_type;
133 :
134 : template <typename T>
135 : struct fill_value_or_no_such_type<T, false> {
136 : using type = NoSuchType;
137 : };
138 :
139 : template <typename T>
140 : struct DoubleWrapper {
141 : static constexpr double value = T::fill_invalid_points_with;
142 : };
143 :
144 : template <typename T>
145 : struct fill_value_or_no_such_type<T, true> {
146 : using type = DoubleWrapper<T>;
147 : };
148 :
149 : template <typename Callback>
150 : struct get_fill_value {
151 : using type = typename fill_value_or_no_such_type<
152 : Callback, has_fill_invalid_points_with_v<Callback>>::type;
153 : };
154 :
155 : template <typename Callbacks>
156 : struct get_invalid_fill {
157 : using type =
158 : tmpl::filter<tmpl::transform<Callbacks, get_fill_value<tmpl::_1>>,
159 : tt::is_a<DoubleWrapper, tmpl::_1>>;
160 : };
161 :
162 : template <typename InterpolationTargetTag, typename TemporalId, typename DbTags>
163 : void fill_invalid_points(const gsl::not_null<db::DataBox<DbTags>*> box,
164 : const TemporalId& temporal_id) {
165 : using callbacks =
166 : typename InterpolationTargetTag::post_interpolation_callbacks;
167 : using invalid_fill = typename get_invalid_fill<callbacks>::type;
168 : if constexpr (tmpl::size<invalid_fill>::value > 0) {
169 : const auto& invalid_indices =
170 : db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(*box);
171 : if (invalid_indices.find(temporal_id) != invalid_indices.end() and
172 : not invalid_indices.at(temporal_id).empty()) {
173 : db::mutate<Tags::IndicesOfInvalidInterpPoints<TemporalId>,
174 : Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
175 : [&temporal_id](
176 : const gsl::not_null<
177 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
178 : indices_of_invalid_points,
179 : const gsl::not_null<std::unordered_map<
180 : TemporalId, Variables<typename InterpolationTargetTag::
181 : vars_to_interpolate_to_target>>*>
182 : vars_dest_all_times) {
183 : auto& vars_dest = vars_dest_all_times->at(temporal_id);
184 : const size_t npts_dest = vars_dest.number_of_grid_points();
185 : const size_t nvars = vars_dest.number_of_independent_components;
186 : for (auto index : indices_of_invalid_points->at(temporal_id)) {
187 : for (size_t v = 0; v < nvars; ++v) {
188 : // clang-tidy: no pointer arithmetic
189 : vars_dest.data()[index + v * npts_dest] = // NOLINT
190 : tmpl::front<invalid_fill>::value;
191 : }
192 : }
193 : // Further functions may test if there are invalid points.
194 : // Clear the invalid points now, since we have filled them.
195 : indices_of_invalid_points->erase(temporal_id);
196 : },
197 : box);
198 : }
199 : } else {
200 : (void)box;
201 : (void)temporal_id;
202 : }
203 : }
204 :
205 : /// Wraps calling the callback function on an InterpolationTarget.
206 : ///
207 : /// First prepares for the callback, then calls the callback,
208 : /// and returns true if the callback is done and
209 : /// false if the callback is not done (e.g. if the callback is an
210 : /// apparent horizon search and it needs another iteration).
211 : ///
212 : /// call_callbacks is called by an Action of InterpolationTarget.
213 : ///
214 : /// Currently two Actions call call_callbacks:
215 : /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
216 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
217 : template <typename InterpolationTargetTag, typename DbTags,
218 : typename Metavariables, typename TemporalId>
219 : bool call_callbacks(
220 : const gsl::not_null<db::DataBox<DbTags>*> box,
221 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
222 : const TemporalId& temporal_id) {
223 : // Before doing anything else, deal with the possibility that some
224 : // of the points might be outside of the Domain.
225 : fill_invalid_points<InterpolationTargetTag>(box, temporal_id);
226 :
227 : // Fill ::Tags::Variables<typename
228 : // InterpolationTargetTag::vars_to_interpolate_to_target>
229 : // with variables from correct temporal_id.
230 : db::mutate_apply<
231 : tmpl::list<::Tags::Variables<
232 : typename InterpolationTargetTag::vars_to_interpolate_to_target>>,
233 : tmpl::list<Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>>(
234 : [&temporal_id](
235 : const gsl::not_null<Variables<
236 : typename InterpolationTargetTag::vars_to_interpolate_to_target>*>
237 : vars,
238 : const std::unordered_map<
239 : TemporalId, Variables<typename InterpolationTargetTag::
240 : vars_to_interpolate_to_target>>&
241 : vars_at_all_times) { *vars = vars_at_all_times.at(temporal_id); },
242 : box);
243 :
244 : // apply_callbacks should return true if we are done with this
245 : // temporal_id. It should return false only if the callback
246 : // calls another `intrp::Action` that needs the volume data at this
247 : // same temporal_id.
248 : return apply_callbacks<InterpolationTargetTag>(box, cache, temporal_id);
249 : }
250 :
251 : /// Frees InterpolationTarget's memory associated with the supplied
252 : /// temporal_id.
253 : ///
254 : /// clean_up_interpolation_target is called by an Action of InterpolationTarget.
255 : ///
256 : /// Currently two Actions call clean_up_interpolation_target:
257 : /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
258 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
259 : template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
260 : void clean_up_interpolation_target(
261 : const gsl::not_null<db::DataBox<DbTags>*> box,
262 : const TemporalId& temporal_id) {
263 : // We are now done with this temporal_id, so we can pop it and
264 : // clean up data associated with it.
265 : db::mutate<Tags::TemporalIds<TemporalId>,
266 : Tags::CompletedTemporalIds<TemporalId>,
267 : Tags::IndicesOfFilledInterpPoints<TemporalId>,
268 : Tags::IndicesOfInvalidInterpPoints<TemporalId>,
269 : Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
270 : [&temporal_id](
271 : const gsl::not_null<std::deque<TemporalId>*> ids,
272 : const gsl::not_null<std::deque<TemporalId>*> completed_ids,
273 : const gsl::not_null<
274 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
275 : indices_of_filled,
276 : const gsl::not_null<
277 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
278 : indices_of_invalid,
279 : const gsl::not_null<std::unordered_map<
280 : TemporalId, Variables<typename InterpolationTargetTag::
281 : vars_to_interpolate_to_target>>*>
282 : interpolated_vars) {
283 : completed_ids->push_back(temporal_id);
284 : ASSERT(std::find(ids->begin(), ids->end(), temporal_id) != ids->end(),
285 : "Temporal id " << temporal_id << " does not exist in ids");
286 : ids->erase(std::remove(ids->begin(), ids->end(), temporal_id),
287 : ids->end());
288 : // We want to keep track of all completed temporal_ids to deal with
289 : // the possibility of late calls to
290 : // AddTemporalIdsToInterpolationTarget. We could keep all
291 : // completed_ids forever, but we probably don't want it to get too
292 : // large, so we limit its size. We assume that
293 : // asynchronous calls to AddTemporalIdsToInterpolationTarget do not span
294 : // more than 1000 temporal_ids.
295 : if (completed_ids->size() > 1000) {
296 : completed_ids->pop_front();
297 : }
298 : indices_of_filled->erase(temporal_id);
299 : indices_of_invalid->erase(temporal_id);
300 : interpolated_vars->erase(temporal_id);
301 : },
302 : box);
303 : }
304 :
305 : /// Returns true if this InterpolationTarget has received data
306 : /// at all its points.
307 : ///
308 : /// have_data_at_all_points is called by an Action of InterpolationTarget.
309 : ///
310 : /// Currently two Actions call have_data_at_all_points:
311 : /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
312 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
313 : template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
314 : bool have_data_at_all_points(const db::DataBox<DbTags>& box,
315 : const TemporalId& temporal_id) {
316 : const auto& filled_indices =
317 : db::get<Tags::IndicesOfFilledInterpPoints<TemporalId>>(box);
318 :
319 : const size_t filled_size = filled_indices.count(temporal_id) > 0
320 : ? filled_indices.at(temporal_id).size()
321 : : 0;
322 : const size_t invalid_size = [&box, &temporal_id]() {
323 : const auto& invalid_indices =
324 : db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(box);
325 : if (invalid_indices.count(temporal_id) > 0) {
326 : return invalid_indices.at(temporal_id).size();
327 : }
328 : return 0_st;
329 : }();
330 : const size_t interp_size =
331 : db::get<Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(box)
332 : .at(temporal_id)
333 : .number_of_grid_points();
334 : return (invalid_size + filled_size == interp_size);
335 : }
336 :
337 : /// Tells an InterpolationTarget that it should interpolate at
338 : /// the supplied temporal_ids. Changes the InterpolationTarget's DataBox
339 : /// accordingly.
340 : ///
341 : /// Returns the temporal_ids that have actually been newly flagged
342 : /// (since some of them may have been flagged already).
343 : ///
344 : /// flag_temporal_ids_for_interpolation is called by an Action
345 : /// of InterpolationTarget
346 : ///
347 : /// Currently one Action calls flag_temporal_ids_for_interpolation:
348 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
349 : template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
350 : std::vector<TemporalId> flag_temporal_ids_for_interpolation(
351 : const gsl::not_null<db::DataBox<DbTags>*> box,
352 : const std::vector<TemporalId>& temporal_ids) {
353 : // We allow this function to be called multiple times with the same
354 : // temporal_ids (e.g. from each element, or from each node of a
355 : // NodeGroup ParallelComponent such as Interpolator). If multiple
356 : // calls occur, we care only about the first one, and ignore the
357 : // others. The first call will often begin interpolation. So if
358 : // multiple calls occur, it is possible that some of them may arrive
359 : // late, even after interpolation has been completed on one or more
360 : // of the temporal_ids (and after that id has already been removed
361 : // from `ids`). If this happens, we don't want to add the
362 : // temporal_ids again. For that reason we keep track of the
363 : // temporal_ids that we have already completed interpolation on. So
364 : // here we do not add any temporal_ids that are already present in
365 : // `ids` or `completed_ids`.
366 : std::vector<TemporalId> new_temporal_ids{};
367 :
368 : db::mutate_apply<tmpl::list<Tags::TemporalIds<TemporalId>>,
369 : tmpl::list<Tags::CompletedTemporalIds<TemporalId>>>(
370 : [&temporal_ids, &new_temporal_ids](
371 : const gsl::not_null<std::deque<TemporalId>*> ids,
372 : const std::deque<TemporalId>& completed_ids) {
373 : for (auto& id : temporal_ids) {
374 : if (std::find(completed_ids.begin(), completed_ids.end(), id) ==
375 : completed_ids.end() and
376 : std::find(ids->begin(), ids->end(), id) == ids->end()) {
377 : ids->push_back(id);
378 : new_temporal_ids.push_back(id);
379 : }
380 : }
381 : },
382 : box);
383 :
384 : return new_temporal_ids;
385 : }
386 :
387 : /// Tells an InterpolationTarget that it should interpolate at
388 : /// the supplied temporal_ids. Changes the InterpolationTarget's DataBox
389 : /// accordingly.
390 : ///
391 : /// Returns the temporal_ids that have actually been newly flagged
392 : /// (since some of them may have been flagged already).
393 : ///
394 : /// flag_temporal_ids_as_pending is called by an Action
395 : /// of InterpolationTarget
396 : ///
397 : /// Currently one Action calls flag_temporal_ids_as_pending:
398 : /// - AddTemporalIdsToInterpolationTarget (called by Events::Interpolate)
399 : template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
400 : std::vector<TemporalId> flag_temporal_ids_as_pending(
401 : const gsl::not_null<db::DataBox<DbTags>*> box,
402 : const std::vector<TemporalId>& temporal_ids) {
403 : // We allow this function to be called multiple times with the same
404 : // temporal_ids (e.g. from each element, or from each node of a
405 : // NodeGroup ParallelComponent such as Interpolator). If multiple
406 : // calls occur, we care only about the first one, and ignore the
407 : // others. The first call will often begin interpolation. So if
408 : // multiple calls occur, it is possible that some of them may arrive
409 : // late, even after interpolation has been completed on one or more
410 : // of the temporal_ids (and after that id has already been removed
411 : // from `ids`). If this happens, we don't want to add the
412 : // temporal_ids again. For that reason we keep track of the
413 : // temporal_ids that we have already completed interpolation on. So
414 : // here we do not add any temporal_ids that are already present in
415 : // `ids` or `completed_ids`.
416 : std::vector<TemporalId> new_temporal_ids{};
417 :
418 : db::mutate_apply<tmpl::list<Tags::PendingTemporalIds<TemporalId>>,
419 : tmpl::list<Tags::TemporalIds<TemporalId>,
420 : Tags::CompletedTemporalIds<TemporalId>>>(
421 : [&temporal_ids, &new_temporal_ids](
422 : const gsl::not_null<std::deque<TemporalId>*> pending_ids,
423 : const std::deque<TemporalId>& ids,
424 : const std::deque<TemporalId>& completed_ids) {
425 : for (auto& id : temporal_ids) {
426 : if (std::find(completed_ids.begin(), completed_ids.end(), id) ==
427 : completed_ids.end() and
428 : std::find(ids.begin(), ids.end(), id) == ids.end() and
429 : std::find(pending_ids->begin(), pending_ids->end(), id) ==
430 : pending_ids->end()) {
431 : pending_ids->push_back(id);
432 : new_temporal_ids.push_back(id);
433 : }
434 : }
435 : },
436 : box);
437 :
438 : return new_temporal_ids;
439 : }
440 :
441 : /// Adds the supplied interpolated variables and offsets to the
442 : /// InterpolationTarget's internal DataBox.
443 : ///
444 : /// Note that the template argument to Variables in vars_src is called
445 : /// InterpolationTargetTag::vars_to_interpolate_to_target. This is a list
446 : /// of tags, and is used for both the interpolated variables (in
447 : /// this function add_received_variables) and for the source variables
448 : /// (in other functions). The source and interpolated quantities are
449 : /// the same set of variables (but at different points).
450 : ///
451 : /// add_received_variables is called by an Action of InterpolationTarget.
452 : ///
453 : /// Currently two Actions call add_received_variables:
454 : /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
455 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
456 : template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
457 : void add_received_variables(
458 : const gsl::not_null<db::DataBox<DbTags>*> box,
459 : const std::vector<Variables<
460 : typename InterpolationTargetTag::vars_to_interpolate_to_target>>&
461 : vars_src,
462 : const std::vector<std::vector<size_t>>& global_offsets,
463 : const TemporalId& temporal_id) {
464 : db::mutate<Tags::IndicesOfFilledInterpPoints<TemporalId>,
465 : Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
466 : [&temporal_id, &vars_src, &global_offsets](
467 : const gsl::not_null<
468 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
469 : indices_of_filled,
470 : const gsl::not_null<std::unordered_map<
471 : TemporalId, Variables<typename InterpolationTargetTag::
472 : vars_to_interpolate_to_target>>*>
473 : vars_dest_all_times) {
474 : auto& vars_dest = vars_dest_all_times->at(temporal_id);
475 : // Here we assume that vars_dest has been allocated to the correct
476 : // size (but could contain garbage, since below we are filling it).
477 : const size_t npts_dest = vars_dest.number_of_grid_points();
478 : const size_t nvars = vars_dest.number_of_independent_components;
479 : for (size_t j = 0; j < global_offsets.size(); ++j) {
480 : const size_t npts_src = global_offsets[j].size();
481 : for (size_t i = 0; i < npts_src; ++i) {
482 : // If a point is on the boundary of two (or more)
483 : // elements, it is possible that we have received data
484 : // for this point from more than one Interpolator.
485 : // This will rarely occur, but it does occur, e.g. when
486 : // a point is exactly on some symmetry
487 : // boundary (such as the x-y plane) and this symmetry
488 : // boundary is exactly the boundary between two
489 : // elements. If this happens, we accept the first
490 : // duplicated point, and we ignore subsequent
491 : // duplicated points. The points are easy to keep track
492 : // of because global_offsets uniquely identifies them.
493 : if ((*indices_of_filled)[temporal_id]
494 : .insert(global_offsets[j][i])
495 : .second) {
496 : for (size_t v = 0; v < nvars; ++v) {
497 : // clang-tidy: no pointer arithmetic
498 : vars_dest.data()[global_offsets[j][i] + // NOLINT
499 : v * npts_dest] = // NOLINT
500 : vars_src[j].data()[i + v * npts_src]; // NOLINT
501 : }
502 : }
503 : }
504 : }
505 : },
506 : box);
507 : }
508 :
509 : /// Computes the block logical coordinates of an InterpolationTarget.
510 : ///
511 : /// block_logical_coords is called by an Action of InterpolationTarget.
512 : ///
513 : /// Currently one Action directly calls this version of block_logical_coords:
514 : /// - InterpolationTargetSendTimeIndepPointsToElements
515 : /// (in InterpolationTarget ActionList)
516 : /// and one Action indirectly calls this version of block_logical_coords:
517 : /// - SendPointsToInterpolator (called by AddTemporalIdsToInterpolationTarget
518 : /// and by FindApparentHorizon)
519 : template <typename InterpolationTargetTag, typename Metavariables,
520 : typename TemporalId>
521 : auto block_logical_coords(
522 : const Parallel::GlobalCache<Metavariables>& cache,
523 : const tnsr::I<
524 : DataVector, Metavariables::volume_dim,
525 : typename InterpolationTargetTag::compute_target_points::frame>& coords,
526 : const TemporalId& temporal_id) {
527 : const auto& domain =
528 : get<domain::Tags::Domain<Metavariables::volume_dim>>(cache);
529 : if constexpr (std::is_same_v<typename InterpolationTargetTag::
530 : compute_target_points::frame,
531 : ::Frame::Grid>) {
532 : // Frame is grid frame, so don't need any FunctionsOfTime,
533 : // whether or not the maps are time_dependent.
534 : return ::block_logical_coordinates(domain, coords);
535 : }
536 :
537 : if (domain.is_time_dependent()) {
538 : if constexpr (Parallel::is_in_global_cache<Metavariables,
539 : domain::Tags::FunctionsOfTime>) {
540 : // Whoever calls block_logical_coords when the maps are
541 : // time-dependent is responsible for ensuring
542 : // that functions_of_time are up to date at temporal_id.
543 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
544 : return ::block_logical_coordinates(
545 : domain, coords,
546 : InterpolationTarget_detail::get_temporal_id_value(temporal_id),
547 : functions_of_time);
548 : } else {
549 : // We error here because the maps are time-dependent, yet
550 : // the cache does not contain FunctionsOfTime. It would be
551 : // nice to make this a compile-time error; however, we want
552 : // the code to compile for the completely time-independent
553 : // case where there are no FunctionsOfTime in the cache at
554 : // all. Unfortunately, checking whether the maps are
555 : // time-dependent is currently not constexpr.
556 : ERROR(
557 : "There is a time-dependent CoordinateMap in at least one "
558 : "of the Blocks, but FunctionsOfTime are not in the "
559 : "GlobalCache. If you intend to use a time-dependent "
560 : "CoordinateMap, please add FunctionsOfTime to the GlobalCache.");
561 : }
562 : }
563 :
564 : // Time-independent case.
565 : return ::block_logical_coordinates(domain, coords);
566 : }
567 :
568 : /// Version of block_logical_coords that computes the interpolation
569 : /// points on the fly. This version is called when the interpolation
570 : /// points are (or might be) time-dependent in the
571 : /// InterpolationTargetTag's frame.
572 : ///
573 : /// This version of block_logical_coordinates is called when there
574 : /// is an Interpolator ParallelComponent.
575 : ///
576 : /// Currently one Action directly calls this version of block_logical_coords:
577 : /// - SendPointsToInterpolator (called by AddTemporalIdsToInterpolationTarget
578 : /// and by FindApparentHorizon)
579 : template <typename InterpolationTargetTag, typename DbTags,
580 : typename Metavariables, typename TemporalId>
581 : auto block_logical_coords(const db::DataBox<DbTags>& box,
582 : const Parallel::GlobalCache<Metavariables>& cache,
583 : const TemporalId& temporal_id) {
584 : return block_logical_coords<InterpolationTargetTag>(
585 : cache,
586 : InterpolationTargetTag::compute_target_points::points(
587 : box, tmpl::type_<Metavariables>{}, temporal_id),
588 : temporal_id);
589 : }
590 :
591 : /// Version of block_logical_coords for when the coords are
592 : /// time-independent.
593 : template <typename InterpolationTargetTag, typename DbTags,
594 : typename Metavariables>
595 : auto block_logical_coords(const db::DataBox<DbTags>& box,
596 : const tmpl::type_<Metavariables>& meta) {
597 : const auto& domain =
598 : db::get<domain::Tags::Domain<Metavariables::volume_dim>>(box);
599 : return ::block_logical_coordinates(
600 : domain, InterpolationTargetTag::compute_target_points::points(box, meta));
601 : }
602 :
603 : /// Initializes InterpolationTarget's variables storage and lists of indices
604 : /// corresponding to the supplied block logical coordinates and `temporal_id`.
605 : ///
606 : /// set_up_interpolation is called by an Action of InterpolationTarget.
607 : ///
608 : /// Currently two Actions call set_up_interpolation:
609 : /// - SendPointsToInterpolator (called by AddTemporalIdsToInterpolationTarget
610 : /// and by FindApparentHorizon)
611 : /// - InterpolationTargetVarsFromElement (called by DgElementArray)
612 : template <typename InterpolationTargetTag, typename DbTags, size_t VolumeDim,
613 : typename TemporalId>
614 : void set_up_interpolation(
615 : const gsl::not_null<db::DataBox<DbTags>*> box,
616 : const TemporalId& temporal_id,
617 : const std::vector<std::optional<
618 : IdPair<domain::BlockId,
619 : tnsr::I<double, VolumeDim, typename ::Frame::BlockLogical>>>>&
620 : block_logical_coords) {
621 : db::mutate<Tags::IndicesOfFilledInterpPoints<TemporalId>,
622 : Tags::IndicesOfInvalidInterpPoints<TemporalId>,
623 : Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
624 : [&block_logical_coords, &temporal_id](
625 : const gsl::not_null<
626 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
627 : indices_of_filled,
628 : const gsl::not_null<
629 : std::unordered_map<TemporalId, std::unordered_set<size_t>>*>
630 : indices_of_invalid_points,
631 : const gsl::not_null<std::unordered_map<
632 : TemporalId, Variables<typename InterpolationTargetTag::
633 : vars_to_interpolate_to_target>>*>
634 : vars_dest_all_times) {
635 : // Because we are sending new points to the interpolator,
636 : // we know that none of these points have been interpolated to,
637 : // so clear the list.
638 : indices_of_filled->erase(temporal_id);
639 :
640 : // Set the indices of invalid points.
641 : indices_of_invalid_points->erase(temporal_id);
642 : for (size_t i = 0; i < block_logical_coords.size(); ++i) {
643 : // The sphere target is optimized specially. Because of this, a
644 : // nullopt in block_logical_coords from the sphere target doesn't
645 : // actually mean the point is invalid. Therefore we ignore this check
646 : // for the sphere target. The downside of this is that we don't catch
647 : // invalid points here.
648 : constexpr bool is_sphere = tt::is_a_v<
649 : TargetPoints::Sphere,
650 : typename InterpolationTargetTag::compute_target_points>;
651 : if (not is_sphere and not block_logical_coords[i].has_value()) {
652 : (*indices_of_invalid_points)[temporal_id].insert(i);
653 : }
654 : }
655 :
656 : // At this point we don't know if vars_dest exists in the map;
657 : // if it doesn't then we want to default construct it.
658 : auto& vars_dest = (*vars_dest_all_times)[temporal_id];
659 :
660 : // We will be filling vars_dest with interpolated data.
661 : // Here we make sure it is allocated to the correct size.
662 : if (vars_dest.number_of_grid_points() != block_logical_coords.size()) {
663 : vars_dest = Variables<
664 : typename InterpolationTargetTag::vars_to_interpolate_to_target>(
665 : block_logical_coords.size());
666 : }
667 : },
668 : box);
669 : }
670 :
671 : CREATE_HAS_TYPE_ALIAS(compute_vars_to_interpolate)
672 : CREATE_HAS_TYPE_ALIAS_V(compute_vars_to_interpolate)
673 :
674 : namespace detail {
675 : template <typename Tag, typename Frame>
676 : using any_index_in_frame_impl =
677 : TensorMetafunctions::any_index_in_frame<typename Tag::type, Frame>;
678 : } // namespace detail
679 :
680 : /// Returns true if any of the tensors in TagList have any of their
681 : /// indices in the given frame.
682 : template <typename TagList, typename Frame>
683 : constexpr bool any_index_in_frame_v =
684 : tmpl::any<TagList, tmpl::bind<detail::any_index_in_frame_impl, tmpl::_1,
685 : Frame>>::value;
686 :
687 : /// Calls compute_vars_to_interpolate to compute
688 : /// InterpolationTargetTag::vars_to_interpolate_to_target from the source
689 : /// variables. Does any frame transformations needed.
690 : template <typename InterpolationTargetTag, typename SourceTags,
691 : typename Metavariables, typename ElementId>
692 : void compute_dest_vars_from_source_vars(
693 : const gsl::not_null<Variables<
694 : typename InterpolationTargetTag::vars_to_interpolate_to_target>*>
695 : dest_vars,
696 : const Variables<SourceTags>& source_vars,
697 : const Domain<Metavariables::volume_dim>& domain,
698 : const Mesh<Metavariables::volume_dim>& mesh, const ElementId& element_id,
699 : const Parallel::GlobalCache<Metavariables>& cache,
700 : const typename InterpolationTargetTag::temporal_id::type& temporal_id) {
701 : const auto& block = domain.blocks().at(element_id.block_id());
702 : if (block.is_time_dependent()) {
703 : // The functions of time are always guaranteed to be
704 : // up-to-date here.
705 : // For interpolation without an Interpolator ParallelComponent,
706 : // this is because the InterpWithoutInterpComponent event will be called
707 : // after the Action that keeps functions of time up to date.
708 : // For interpolation with an Interpolator ParallelComponent,
709 : // this is because the functions of time are made up to date before
710 : // calling SendPointsToInterpolator.
711 : if constexpr (any_index_in_frame_v<SourceTags, Frame::Inertial> and
712 : any_index_in_frame_v<typename InterpolationTargetTag::
713 : vars_to_interpolate_to_target,
714 : Frame::Grid>) {
715 : // Need to do frame transformations to Grid frame.
716 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
717 : ElementMap<3, ::Frame::Grid> map_logical_to_grid{
718 : element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
719 : const auto logical_coords = logical_coordinates(mesh);
720 : const auto time =
721 : InterpolationTarget_detail::get_temporal_id_value(temporal_id);
722 : const auto jac_logical_to_grid =
723 : map_logical_to_grid.jacobian(logical_coords);
724 : const auto invjac_logical_to_grid =
725 : map_logical_to_grid.inv_jacobian(logical_coords);
726 : const auto [inertial_coords, invjac_grid_to_inertial,
727 : jac_grid_to_inertial, inertial_mesh_velocity] =
728 : block.moving_mesh_grid_to_inertial_map()
729 : .coords_frame_velocity_jacobians(
730 : map_logical_to_grid(logical_coords), time, functions_of_time);
731 : InterpolationTargetTag::compute_vars_to_interpolate::apply(
732 : dest_vars, source_vars, mesh, jac_grid_to_inertial,
733 : invjac_grid_to_inertial, jac_logical_to_grid, invjac_logical_to_grid,
734 : inertial_mesh_velocity,
735 : tnsr::I<DataVector, 3, Frame::Grid>{
736 : DataVector{get<0, 0>(jac_logical_to_grid).size(), 0.0}});
737 : } else if constexpr (any_index_in_frame_v<SourceTags, Frame::Inertial> and
738 : any_index_in_frame_v<typename InterpolationTargetTag::
739 : vars_to_interpolate_to_target,
740 : Frame::Distorted>) {
741 : // Need to do frame transformations to Distorted frame.
742 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
743 : ASSERT(block.has_distorted_frame(),
744 : "Cannot interpolate to distorted frame in a block that does not "
745 : "have a distorted frame");
746 : const domain::CoordinateMaps::Composition
747 : element_logical_to_distorted_map{
748 : domain::element_to_block_logical_map(element_id),
749 : block.moving_mesh_logical_to_grid_map().get_clone(),
750 : block.moving_mesh_grid_to_distorted_map().get_clone()};
751 : const domain::CoordinateMaps::Composition element_logical_to_grid_map{
752 : domain::element_to_block_logical_map(element_id),
753 : block.moving_mesh_logical_to_grid_map().get_clone()};
754 : const auto logical_coords = logical_coordinates(mesh);
755 : const auto time =
756 : InterpolationTarget_detail::get_temporal_id_value(temporal_id);
757 : const auto [inertial_coords, invjac_distorted_to_inertial,
758 : jac_distorted_to_inertial,
759 : distorted_to_inertial_mesh_velocity] =
760 : block.moving_mesh_distorted_to_inertial_map()
761 : .coords_frame_velocity_jacobians(
762 : element_logical_to_distorted_map(logical_coords, time,
763 : functions_of_time),
764 : time, functions_of_time);
765 : const auto grid_to_distorted_mesh_velocity =
766 : get<3>(block.moving_mesh_grid_to_distorted_map()
767 : .coords_frame_velocity_jacobians(
768 : element_logical_to_grid_map(logical_coords, time,
769 : functions_of_time),
770 : time, functions_of_time));
771 : InterpolationTargetTag::compute_vars_to_interpolate::apply(
772 : dest_vars, source_vars, mesh, jac_distorted_to_inertial,
773 : invjac_distorted_to_inertial,
774 : element_logical_to_distorted_map.jacobian(logical_coords, time,
775 : functions_of_time),
776 : element_logical_to_distorted_map.inv_jacobian(logical_coords, time,
777 : functions_of_time),
778 : distorted_to_inertial_mesh_velocity, grid_to_distorted_mesh_velocity);
779 : } else {
780 : // No frame transformations needed.
781 : InterpolationTargetTag::compute_vars_to_interpolate::apply(
782 : dest_vars, source_vars, mesh);
783 : }
784 : } else {
785 : // No frame transformations needed, since the maps are time-independent
786 : // and therefore all the frames are the same.
787 : //
788 : // Sometimes dest_vars and source_vars have different Frame tags
789 : // in the time-independent case, even though the frames are really
790 : // the same. The source vars should all be in the Inertial frame,
791 : // so we create a new non-owning Variables called
792 : // dest_vars_in_inertial_frame that points to dest_vars but is
793 : // tagged as the inertial frame, and we pass
794 : // dest_vars_in_inertial_frame to
795 : // compute_vars_to_interpolate::apply.
796 : using dest_vars_tags =
797 : typename InterpolationTargetTag::vars_to_interpolate_to_target;
798 : using dest_tags_in_inertial_frame =
799 : TensorMetafunctions::replace_frame_in_taglist<dest_vars_tags,
800 : ::Frame::Inertial>;
801 : Variables<dest_tags_in_inertial_frame> dest_vars_in_inertial_frame(
802 : dest_vars->data(), dest_vars->size());
803 : InterpolationTargetTag::compute_vars_to_interpolate::apply(
804 : make_not_null(&dest_vars_in_inertial_frame), source_vars, mesh);
805 : }
806 : }
807 :
808 : } // namespace InterpolationTarget_detail
809 : } // namespace intrp
|
