InterpolationTargetDetail.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <memory>
8 #include <unordered_set>
9 #include <utility>
10 #include <vector>
11 
14 #include "DataStructures/IdPair.hpp"
15 #include "DataStructures/VariablesTag.hpp"
16 #include "Domain/BlockLogicalCoordinates.hpp"
17 #include "Domain/Tags.hpp"
19 #include "Parallel/Invoke.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/Literals.hpp"
22 #include "Utilities/Requires.hpp"
23 #include "Utilities/TMPL.hpp"
24 #include "Utilities/TaggedTuple.hpp"
25 #include "Utilities/TypeTraits.hpp"
26 #include "Utilities/TypeTraits/CreateHasStaticMemberVariable.hpp"
27 
28 /// \cond
29 // IWYU pragma: no_forward_declare db::DataBox
30 namespace intrp {
31 template <class Metavariables>
32 struct Interpolator;
33 template <typename Metavariables, typename InterpolationTargetTag>
34 class InterpolationTarget;
35 namespace Actions {
36 template <typename InterpolationTargetTag>
37 struct CleanUpInterpolator;
38 template <typename InterpolationTargetTag>
39 struct ReceivePoints;
40 } // namespace Actions
41 namespace Tags {
42 template <typename TemporalId>
43 struct IndicesOfFilledInterpPoints;
44 template <typename TemporalId>
45 struct IndicesOfInvalidInterpPoints;
46 template <typename InterpolationTargetTag, typename TemporalId>
47 struct InterpolatedVars;
48 template <typename TemporalId>
49 struct CompletedTemporalIds;
50 template <typename TemporalId>
51 struct TemporalIds;
52 } // namespace Tags
53 } // namespace intrp
54 template <typename TagsList>
55 struct Variables;
56 /// \endcond
57 
58 namespace intrp {
59 
60 namespace InterpolationTarget_detail {
61 
62 // apply_callback accomplishes the overload for the
63 // two signatures of callback functions.
64 // Uses SFINAE on return type.
65 template <typename T, typename DbTags, typename Metavariables,
66  typename TemporalId>
67 auto apply_callback(
70  const TemporalId& temporal_id) noexcept
71  -> decltype(T::post_interpolation_callback::apply(box, cache, temporal_id),
72  bool()) {
73  return T::post_interpolation_callback::apply(box, cache, temporal_id);
74 }
75 
76 template <typename T, typename DbTags, typename Metavariables,
77  typename TemporalId>
78 auto apply_callback(
81  const TemporalId& temporal_id) noexcept
82  -> decltype(T::post_interpolation_callback::apply(*box, *cache,
83  temporal_id),
84  bool()) {
85  T::post_interpolation_callback::apply(*box, *cache, temporal_id);
86  // For the simpler callback function, we will always clean up volume data, so
87  // we return true here.
88  return true;
89 }
90 
91 CREATE_HAS_STATIC_MEMBER_VARIABLE(fill_invalid_points_with)
92 CREATE_HAS_STATIC_MEMBER_VARIABLE_V(fill_invalid_points_with)
93 
94 // Fills invalid points with some constant value.
95 template <typename InterpolationTargetTag, typename TemporalId, typename DbTags,
96  Requires<not has_fill_invalid_points_with_v<
97  typename InterpolationTargetTag::post_interpolation_callback>> =
98  nullptr>
99 void fill_invalid_points(const gsl::not_null<db::DataBox<DbTags>*> /*box*/,
100  const TemporalId& /*temporal_id*/) noexcept {}
101 
102 template <typename InterpolationTargetTag, typename TemporalId, typename DbTags,
103  Requires<has_fill_invalid_points_with_v<
104  typename InterpolationTargetTag::post_interpolation_callback>> =
105  nullptr>
106 void fill_invalid_points(const gsl::not_null<db::DataBox<DbTags>*> box,
107  const TemporalId& temporal_id) noexcept {
108  const auto& invalid_indices =
109  db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(*box);
110  if (invalid_indices.find(temporal_id) != invalid_indices.end() and
111  not invalid_indices.at(temporal_id).empty()) {
112  db::mutate<Tags::IndicesOfInvalidInterpPoints<TemporalId>,
113  Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
114  box, [&temporal_id](
116  TemporalId, std::unordered_set<size_t>>*>
117  indices_of_invalid_points,
119  TemporalId, Variables<typename InterpolationTargetTag::
120  vars_to_interpolate_to_target>>*>
121  vars_dest_all_times) noexcept {
122  auto& vars_dest = vars_dest_all_times->at(temporal_id);
123  const size_t npts_dest = vars_dest.number_of_grid_points();
124  const size_t nvars = vars_dest.number_of_independent_components;
125  for (auto index : indices_of_invalid_points->at(temporal_id)) {
126  for (size_t v = 0; v < nvars; ++v) {
127  // clang-tidy: no pointer arithmetic
128  vars_dest.data()[index + v * npts_dest] = // NOLINT
129  InterpolationTargetTag::post_interpolation_callback::
130  fill_invalid_points_with;
131  }
132  }
133  // Further functions may test if there are invalid points.
134  // Clear the invalid points now, since we have filled them.
135  indices_of_invalid_points->erase(temporal_id);
136  });
137  }
138 }
139 
140 /// Wraps calling the callback function on an InterpolationTarget.
141 ///
142 /// First prepares for the callback, then calls the callback,
143 /// and returns true if the callback is done and
144 /// false if the callback is not done (e.g. if the callback is an
145 /// apparent horizon search and it needs another iteration).
146 ///
147 /// call_callback is called by an Action of InterpolationTarget.
148 ///
149 /// Currently two Actions call call_callback:
150 /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
151 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
152 template <typename InterpolationTargetTag, typename DbTags,
153  typename Metavariables, typename TemporalId>
154 bool call_callback(
157  const TemporalId& temporal_id) noexcept {
158  // Before doing anything else, deal with the possibility that some
159  // of the points might be outside of the Domain.
160  fill_invalid_points<InterpolationTargetTag>(box, temporal_id);
161 
162  // Fill ::Tags::Variables<typename
163  // InterpolationTargetTag::vars_to_interpolate_to_target>
164  // with variables from correct temporal_id.
166  tmpl::list<::Tags::Variables<
167  typename InterpolationTargetTag::vars_to_interpolate_to_target>>,
168  tmpl::list<Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>>(
169  [&temporal_id](
170  const gsl::not_null<Variables<
171  typename InterpolationTargetTag::vars_to_interpolate_to_target>*>
172  vars,
173  const std::unordered_map<
174  TemporalId, Variables<typename InterpolationTargetTag::
175  vars_to_interpolate_to_target>>&
176  vars_at_all_times) noexcept {
177  *vars = vars_at_all_times.at(temporal_id);
178  },
179  box);
180 
181  // apply_callback should return true if we are done with this
182  // temporal_id. It should return false only if the callback
183  // calls another `intrp::Action` that needs the volume data at this
184  // same temporal_id.
185  return apply_callback<InterpolationTargetTag>(box, cache, temporal_id);
186 }
187 
188 /// Frees InterpolationTarget's memory associated with the supplied
189 /// temporal_id.
190 ///
191 /// clean_up_interpolation_target is called by an Action of InterpolationTarget.
192 ///
193 /// Currently two Actions call clean_up_interpolation_target:
194 /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
195 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
196 template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
197 void clean_up_interpolation_target(
199  const TemporalId& temporal_id) noexcept {
200  // We are now done with this temporal_id, so we can pop it and
201  // clean up data associated with it.
202  db::mutate<Tags::TemporalIds<TemporalId>,
203  Tags::CompletedTemporalIds<TemporalId>,
204  Tags::IndicesOfFilledInterpPoints<TemporalId>,
205  Tags::IndicesOfInvalidInterpPoints<TemporalId>,
206  Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
207  box, [&temporal_id](
209  const gsl::not_null<std::deque<TemporalId>*> completed_ids,
210  const gsl::not_null<
212  indices_of_filled,
213  const gsl::not_null<
215  indices_of_invalid,
217  TemporalId, Variables<typename InterpolationTargetTag::
218  vars_to_interpolate_to_target>>*>
219  interpolated_vars) noexcept {
220  completed_ids->push_back(temporal_id);
221  ASSERT(std::find(ids->begin(), ids->end(), temporal_id) != ids->end(),
222  "Temporal id " << temporal_id << " does not exist in ids");
223  ids->erase(std::remove(ids->begin(), ids->end(), temporal_id),
224  ids->end());
225  // We want to keep track of all completed temporal_ids to deal with
226  // the possibility of late calls to
227  // AddTemporalIdsToInterpolationTarget. We could keep all
228  // completed_ids forever, but we probably don't want it to get too
229  // large, so we limit its size. We assume that
230  // asynchronous calls to AddTemporalIdsToInterpolationTarget do not span
231  // more than 1000 temporal_ids.
232  if (completed_ids->size() > 1000) {
233  completed_ids->pop_front();
234  }
235  indices_of_filled->erase(temporal_id);
236  indices_of_invalid->erase(temporal_id);
237  interpolated_vars->erase(temporal_id);
238  });
239 }
240 
241 /// Returns true if this InterpolationTarget has received data
242 /// at all its points.
243 ///
244 /// have_data_at_all_points is called by an Action of InterpolationTarget.
245 ///
246 /// Currently two Actions call have_data_at_all_points:
247 /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
248 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
249 template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
250 bool have_data_at_all_points(const db::DataBox<DbTags>& box,
251  const TemporalId& temporal_id) noexcept {
252  const size_t filled_size =
253  db::get<Tags::IndicesOfFilledInterpPoints<TemporalId>>(box)
254  .at(temporal_id)
255  .size();
256  const size_t invalid_size = [&box, &temporal_id ]() noexcept {
257  const auto& invalid_indices =
258  db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(box);
259  if (invalid_indices.count(temporal_id) > 0) {
260  return invalid_indices.at(temporal_id).size();
261  }
262  return 0_st;
263  }
264  ();
265  const size_t interp_size =
266  db::get<Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(box)
267  .at(temporal_id)
268  .number_of_grid_points();
269  return (invalid_size + filled_size == interp_size);
270 }
271 
272 /// Tells an InterpolationTarget that it should interpolate at
273 /// the supplied temporal_ids. Changes the InterpolationTarget's DataBox
274 /// accordingly.
275 ///
276 /// Returns the temporal_ids that have actually been newly flagged
277 /// (since some of them may have been flagged already).
278 ///
279 /// flag_temporal_ids_for_interpolation is called by an Action
280 /// of InterpolationTarget
281 ///
282 /// Currently two Actions call flag_temporal_ids_for_interpolation:
283 /// - AddTemporalIdsToInterpolationTarget (called by Events::Interpolate)
284 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
285 template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
286 std::vector<TemporalId> flag_temporal_ids_for_interpolation(
288  const std::vector<TemporalId>& temporal_ids) noexcept {
289  // We allow this function to be called multiple times with the same
290  // temporal_ids (e.g. from each element, or from each node of a
291  // NodeGroup ParallelComponent such as Interpolator). If multiple
292  // calls occur, we care only about the first one, and ignore the
293  // others. The first call will often begin interpolation. So if
294  // multiple calls occur, it is possible that some of them may arrive
295  // late, even after interpolation has been completed on one or more
296  // of the temporal_ids (and after that id has already been removed
297  // from `ids`). If this happens, we don't want to add the
298  // temporal_ids again. For that reason we keep track of the
299  // temporal_ids that we have already completed interpolation on. So
300  // here we do not add any temporal_ids that are already present in
301  // `ids` or `completed_ids`.
302  std::vector<TemporalId> new_temporal_ids{};
303 
304  db::mutate_apply<tmpl::list<Tags::TemporalIds<TemporalId>>,
305  tmpl::list<Tags::CompletedTemporalIds<TemporalId>>>(
306  [&temporal_ids, &new_temporal_ids ](
308  const std::deque<TemporalId>& completed_ids) noexcept {
309  for (auto& id : temporal_ids) {
310  if (std::find(completed_ids.begin(), completed_ids.end(), id) ==
311  completed_ids.end() and
312  std::find(ids->begin(), ids->end(), id) == ids->end()) {
313  ids->push_back(id);
314  new_temporal_ids.push_back(id);
315  }
316  }
317  },
318  box);
319 
320  return new_temporal_ids;
321 }
322 
323 /// Adds the supplied interpolated variables and offsets to the
324 /// InterpolationTarget's internal DataBox.
325 ///
326 /// Note that the template argument to Variables in vars_src is called
327 /// InterpolationTargetTag::vars_to_interpolate_to_target. This is a list
328 /// of tags, and is used for both the interpolated variables (in
329 /// this function add_received_variables) and for the source variables
330 /// (in other functions). The source and interpolated quantities are
331 /// the same set of variables (but at different points).
332 ///
333 /// add_received_variables is called by an Action of InterpolationTarget.
334 ///
335 /// Currently two Actions call add_received_variables:
336 /// - InterpolationTargetReceiveVars (called by Interpolator ParallelComponent)
337 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
338 template <typename InterpolationTargetTag, typename DbTags, typename TemporalId>
339 void add_received_variables(
341  const std::vector<Variables<
342  typename InterpolationTargetTag::vars_to_interpolate_to_target>>&
343  vars_src,
344  const std::vector<std::vector<size_t>>& global_offsets,
345  const TemporalId& temporal_id) noexcept {
346  db::mutate<Tags::IndicesOfFilledInterpPoints<TemporalId>,
347  Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
348  box, [&temporal_id, &vars_src, &global_offsets ](
349  const gsl::not_null<
351  indices_of_filled,
353  TemporalId, Variables<typename InterpolationTargetTag::
354  vars_to_interpolate_to_target>>*>
355  vars_dest_all_times) noexcept {
356  auto& vars_dest = vars_dest_all_times->at(temporal_id);
357  // Here we assume that vars_dest has been allocated to the correct
358  // size (but could contain garbage, since below we are filling it).
359  const size_t npts_dest = vars_dest.number_of_grid_points();
360  const size_t nvars = vars_dest.number_of_independent_components;
361  for (size_t j = 0; j < global_offsets.size(); ++j) {
362  const size_t npts_src = global_offsets[j].size();
363  for (size_t i = 0; i < npts_src; ++i) {
364  // If a point is on the boundary of two (or more)
365  // elements, it is possible that we have received data
366  // for this point from more than one Interpolator.
367  // This will rarely occur, but it does occur, e.g. when
368  // a point is exactly on some symmetry
369  // boundary (such as the x-y plane) and this symmetry
370  // boundary is exactly the boundary between two
371  // elements. If this happens, we accept the first
372  // duplicated point, and we ignore subsequent
373  // duplicated points. The points are easy to keep track
374  // of because global_offsets uniquely identifies them.
375  if ((*indices_of_filled)[temporal_id]
376  .insert(global_offsets[j][i])
377  .second) {
378  for (size_t v = 0; v < nvars; ++v) {
379  // clang-tidy: no pointer arithmetic
380  vars_dest.data()[global_offsets[j][i] + // NOLINT
381  v * npts_dest] = // NOLINT
382  vars_src[j].data()[i + v * npts_src]; // NOLINT
383  }
384  }
385  }
386  }
387  });
388 }
389 
390 /// Computes the block logical coordinates of an InterpolationTarget.
391 ///
392 /// block_logical_coords is called by an Action of InterpolationTarget.
393 ///
394 /// Currently two Actions call block_logical_coords:
395 /// - SendPointsToInterpolator (called by AddTemporalIdsToInterpolationTarget
396 /// and by FindApparentHorizon)
397 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
398 /// - InterpolationTargetSendTimeIndepPointsToElements
399 /// (in InterpolationTarget ActionList)
400 template <typename InterpolationTargetTag, typename DbTags,
401  typename Metavariables, typename TemporalId>
402 auto block_logical_coords(const db::DataBox<DbTags>& box,
403  const tmpl::type_<Metavariables>& meta,
404  const TemporalId& temporal_id) noexcept {
405  const auto& domain =
406  db::get<domain::Tags::Domain<Metavariables::volume_dim>>(box);
408  domain, InterpolationTargetTag::compute_target_points::points(
409  box, meta, temporal_id));
410 }
411 
412 /// This is a version of block_logical_coords for when the coords
413 /// are time-independent.
414 template <typename InterpolationTargetTag, typename DbTags,
415  typename Metavariables>
416 auto block_logical_coords(const db::DataBox<DbTags>& box,
417  const tmpl::type_<Metavariables>& meta) noexcept {
418  const auto& domain =
419  db::get<domain::Tags::Domain<Metavariables::volume_dim>>(box);
421  domain, InterpolationTargetTag::compute_target_points::points(box, meta));
422 }
423 
424 /// Initializes InterpolationTarget's variables storage and lists of indices
425 /// corresponding to the supplied block logical coordinates and `temporal_id`.
426 ///
427 /// set_up_interpolation is called by an Action of InterpolationTarget.
428 ///
429 /// Currently two Actions call set_up_interpolation:
430 /// - SendPointsToInterpolator (called by AddTemporalIdsToInterpolationTarget
431 /// and by FindApparentHorizon)
432 /// - InterpolationTargetVarsFromElement (called by DgElementArray)
433 template <typename InterpolationTargetTag, typename DbTags, size_t VolumeDim,
434  typename TemporalId>
435 void set_up_interpolation(
437  const TemporalId& temporal_id,
438  const std::vector<boost::optional<
440  tnsr::I<double, VolumeDim, typename ::Frame::Logical>>>>&
441  block_logical_coords) noexcept {
442  db::mutate<Tags::IndicesOfFilledInterpPoints<TemporalId>,
443  Tags::IndicesOfInvalidInterpPoints<TemporalId>,
444  Tags::InterpolatedVars<InterpolationTargetTag, TemporalId>>(
445  box, [&block_logical_coords, &temporal_id ](
446  const gsl::not_null<
448  indices_of_filled,
449  const gsl::not_null<
451  indices_of_invalid_points,
453  TemporalId, Variables<typename InterpolationTargetTag::
454  vars_to_interpolate_to_target>>*>
455  vars_dest_all_times) noexcept {
456  // Because we are sending new points to the interpolator,
457  // we know that none of these points have been interpolated to,
458  // so clear the list.
459  indices_of_filled->erase(temporal_id);
460 
461  // Set the indices of invalid points.
462  indices_of_invalid_points->erase(temporal_id);
463  for (size_t i = 0; i < block_logical_coords.size(); ++i) {
464  if (not block_logical_coords[i]) {
465  (*indices_of_invalid_points)[temporal_id].insert(i);
466  }
467  }
468 
469  // At this point we don't know if vars_dest exists in the map;
470  // if it doesn't then we want to default construct it.
471  auto& vars_dest = (*vars_dest_all_times)[temporal_id];
472 
473  // We will be filling vars_dest with interpolated data.
474  // Here we make sure it is allocated to the correct size.
475  if (vars_dest.number_of_grid_points() != block_logical_coords.size()) {
476  vars_dest = Variables<
477  typename InterpolationTargetTag::vars_to_interpolate_to_target>(
478  block_logical_coords.size());
479  }
480  });
481 }
482 
483 } // namespace InterpolationTarget_detail
484 } // namespace intrp
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:124
DataBoxTag.hpp
utility
unordered_set
Literals.hpp
Tags.hpp
vector
Tags::Variables
Definition: VariablesTag.hpp:21
IdPair
A data structure that contains an ID and data associated with that ID.
Definition: IdPair.hpp:16
CREATE_HAS_STATIC_MEMBER_VARIABLE
#define CREATE_HAS_STATIC_MEMBER_VARIABLE(CONSTEXPR_NAME)
Generate a type trait to check if a class has a static constexpr variable, optionally also checking i...
Definition: CreateHasStaticMemberVariable.hpp:28
db::apply
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1443
block_logical_coordinates
std::vector< block_logical_coord_holder< Dim > > block_logical_coordinates(const Domain< Dim > &domain, const tnsr::I< DataVector, Dim, Frame::Inertial > &x) noexcept
Definition: BlockLogicalCoordinates.cpp:25
DataBox.hpp
db::mutate_apply
constexpr auto mutate_apply(F &&f, const gsl::not_null< DataBox< BoxTags > * > box, Args &&... args) noexcept(DataBox_detail::check_mutate_apply_mutate_tags(BoxTags{}, MutateTags{}) and DataBox_detail::check_mutate_apply_argument_tags(BoxTags{}, ArgumentTags{}) and noexcept(DataBox_detail::mutate_apply(f, box, MutateTags{}, ArgumentTags{}, std::forward< Args >(args)...)))
Apply the invokable f mutating items MutateTags and taking as additional arguments ArgumentTags and a...
Definition: DataBox.hpp:1659
cstddef
std::deque
memory
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Gsl.hpp
Requires.hpp
domain::BlockId
Index a block of the computational domain.
Definition: BlockId.hpp:21
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
std::unordered_map
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
TMPL.hpp
ConstGlobalCache.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182