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