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