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