ReconstructionCommunication.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <array>
8 #include <boost/functional/hash.hpp>
9 #include <cstddef>
10 #include <iterator>
11 #include <limits>
12 #include <map>
13 #include <tuple>
14 #include <utility>
15 #include <vector>
16 
19 #include "DataStructures/DataVector.hpp"
20 #include "DataStructures/FixedHashMap.hpp"
21 #include "DataStructures/Index.hpp"
24 #include "DataStructures/VariablesTag.hpp"
28 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
29 #include "Domain/Structure/OrientationMapHelpers.hpp"
30 #include "Domain/Tags.hpp"
31 #include "Evolution/DgSubcell/ActiveGrid.hpp"
32 #include "Evolution/DgSubcell/NeighborData.hpp"
33 #include "Evolution/DgSubcell/Projection.hpp"
34 #include "Evolution/DgSubcell/RdmpTci.hpp"
35 #include "Evolution/DgSubcell/SliceData.hpp"
36 #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
37 #include "Evolution/DgSubcell/Tags/Inactive.hpp"
38 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
39 #include "Evolution/DgSubcell/Tags/NeighborData.hpp"
40 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
41 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
42 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
45 #include "Parallel/AlgorithmMetafunctions.hpp"
46 #include "Parallel/GlobalCache.hpp"
47 #include "Time/Tags.hpp"
48 #include "Time/TimeStepId.hpp"
50 #include "Utilities/Gsl.hpp"
51 #include "Utilities/MakeArray.hpp"
52 #include "Utilities/TMPL.hpp"
53 
55 /*!
56  * \brief Sets the local data from the relaxed discrete maximum principle
57  * troubled-cell indicator and sends ghost zone data to neighboring elements.
58  *
59  * The action proceeds as follows:
60  *
61  * 1. Computes the maximum and minimum of each evolved variable, which is used
62  * by the relaxed discrete maximum principle troubled-cell indicator.
63  * 2. Determine in which directions we have neighbors
64  * 3. Slice the _conserved_ variables to send to our neighbors for ghost zones
65  * 4. Send the ghost zone data, appending the max/min for the TCI at the end of
66  * the `std::vector<double>` we are sending.
67  *
68  * Some notes:
69  * - In the future we will need to experiment with sending (some) primitive
70  * variables to the neighbors. This will end up depending on the
71  * reconstruction scheme used, so it'll take a bit of infrastructure
72  * development.
73  * - In the future we will need to send the cell-centered fluxes to do
74  * high-order FD without additional reconstruction being necessary.
75  *
76  * GlobalCache:
77  * - Uses:
78  * - `ParallelComponent` proxy
79  *
80  * DataBox:
81  * - Uses:
82  * - `domain::Tags::Mesh<Dim>`
83  * - `subcell::Tags::Mesh<Dim>`
84  * - `domain::Tags::Element<Dim>`
85  * - `Tags::TimeStepId`
86  * - `Tags::Next<Tags::TimeStepId>`
87  * - `subcell::Tags::ActiveGrid`
88  * - `System::variables_tag`
89  * - Adds: nothing
90  * - Removes: nothing
91  * - Modifies:
92  * - `subcell::Tags::NeighborDataForReconstructionAndRdmpTci<Dim>`
93  */
94 template <size_t Dim, typename GhostDataMutator>
96  using inbox_tags = tmpl::list<
98 
99  template <typename DbTags, typename... InboxTags, typename ArrayIndex,
100  typename ActionList, typename ParallelComponent,
101  typename Metavariables>
102  static std::tuple<db::DataBox<DbTags>&&> apply(
103  db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
105  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
106  const ParallelComponent* const /*meta*/) noexcept {
107  using variables_tag = typename Metavariables::system::variables_tag;
108 
109  ASSERT(db::get<Tags::ActiveGrid>(box) == ActiveGrid::Subcell,
110  "The SendDataForReconstruction action can only be called when "
111  "Subcell is the active scheme.");
112 
113  db::mutate<Tags::NeighborDataForReconstructionAndRdmpTci<Dim>>(
114  make_not_null(&box),
115  [](const auto neighbor_data_ptr, const auto& active_vars) noexcept {
116  auto [max_of_vars, min_of_vars] =
117  rdmp_max_min(active_vars, {}, false);
118 
119  // Clear the previous neighbor data and add current local data
120  neighbor_data_ptr->clear();
121  (*neighbor_data_ptr)[std::pair{
124  NeighborData{{}, std::move(max_of_vars), std::move(min_of_vars)};
125  },
126  db::get<variables_tag>(box));
127 
128  const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
129  const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
130  const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
131  const size_t ghost_zone_size =
132  Metavariables::SubcellOptions::ghost_zone_size(box);
133  // Note: we currently always send the conserved variables, but will need to
134  // generalize that to allow the user to choose what to send. For example,
135  // with relativistic hydrodynamics it's desirable to reconstruct
136  // u_i = W v_i = \alpha u^0 v_i
137  // since this guarantees that the Lorentz factor is physical after
138  // reconstruction.
139  DirectionMap<Dim, bool> directions_to_slice{};
140  for (const auto& direction_neighbors : element.neighbors()) {
141  if (direction_neighbors.second.size() == 0) {
142  directions_to_slice[direction_neighbors.first] = false;
143  } else {
144  directions_to_slice[direction_neighbors.first] = true;
145  }
146  }
147  // Optimization note: could save a copy+allocation if we moved
148  // all_sliced_data when possible before sending.
149  const DirectionMap<Dim, std::vector<double>> all_sliced_data = slice_data(
150  db::mutate_apply(GhostDataMutator{}, make_not_null(&box)),
151  subcell_mesh.extents(), ghost_zone_size, directions_to_slice);
152 
153  auto& receiver_proxy =
154  Parallel::get_parallel_component<ParallelComponent>(cache);
155  const NeighborData& local_neighbor_data =
156  db::get<Tags::NeighborDataForReconstructionAndRdmpTci<Dim>>(box).at(
159  const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
160  const TimeStepId& next_time_step_id = [&box]() noexcept {
161  if (Metavariables::local_time_stepping) {
162  return db::get<::Tags::Next<::Tags::TimeStepId>>(box);
163  } else {
164  return db::get<::Tags::TimeStepId>(box);
165  }
166  }();
167 
168  // Compute and send actual variables
169  for (const auto& [direction, neighbors_in_direction] :
170  element.neighbors()) {
171  const auto& orientation = neighbors_in_direction.orientation();
172  const auto direction_from_neighbor = orientation(direction.opposite());
173  // If we are periodic in this direction don't send any data.
174  ASSERT(neighbors_in_direction.size() == 1,
175  "AMR is not yet supported when using DG-subcell. Note that this "
176  "condition could be relaxed to support AMR only where the "
177  "evolution is using DG without any changes to subcell.");
178 
179  for (const ElementId<Dim>& neighbor : neighbors_in_direction) {
180  std::vector<double> subcell_data_to_send{};
181  if (not orientation.is_aligned()) {
182  std::array<size_t, Dim> slice_extents{};
183  for (size_t d = 0; d < Dim; ++d) {
184  gsl::at(slice_extents, d) = subcell_mesh.extents(d);
185  }
186  gsl::at(slice_extents, direction.dimension()) = ghost_zone_size;
187 
188  subcell_data_to_send =
189  orient_variables(all_sliced_data.at(direction),
190  Index<Dim>{slice_extents}, orientation);
191  } else {
192  subcell_data_to_send = all_sliced_data.at(direction);
193  }
194  subcell_data_to_send.insert(
195  subcell_data_to_send.end(),
196  local_neighbor_data.max_variables_values.cbegin(),
197  local_neighbor_data.max_variables_values.cend());
198  subcell_data_to_send.insert(
199  subcell_data_to_send.end(),
200  local_neighbor_data.min_variables_values.cbegin(),
201  local_neighbor_data.min_variables_values.cend());
202 
205  data{dg_mesh.slice_away(direction.dimension()),
206  std::move(subcell_data_to_send), std::nullopt,
207  next_time_step_id};
208 
211  receiver_proxy[neighbor], time_step_id,
212  std::pair{std::pair{direction_from_neighbor, element.id()},
213  std::move(data)});
214  }
215  }
216  return {std::move(box)};
217  }
218 };
219 
220 /*!
221  * \brief Receive the subcell data from our neighbor, and accumulate the data
222  * from the relaxed discrete maximum principle troubled-cell indicator.
223  *
224  * Note:
225  * - Since we only care about the min/max over all neighbors and ourself at the
226  * past time, we accumulate all data immediately into the self `NeighborData`.
227  * - If the neighbor is using DG and therefore sends boundary correction data
228  * then that is added into the `evolution::dg::Tags::MortarData` tag
229  * - The next `TimeStepId` is recorded, but we do not yet support local time
230  * stepping.
231  * - This action will never care about what variables are sent for
232  * reconstruction. It is only responsible for receiving the data and storing
233  * it in the `NeighborData`.
234  *
235  * GlobalCache:
236  * -Uses: nothing
237  *
238  * DataBox:
239  * - Uses:
240  * - `domain::Tags::Element<Dim>`
241  * - `Tags::TimeStepId`
242  * - `domain::Tags::Mesh<Dim>`
243  * - `subcell::Tags::Mesh<Dim>`
244  * - `domain::Tags::Element<Dim>`
245  * - `Tags::Next<Tags::TimeStepId>`
246  * - `subcell::Tags::ActiveGrid`
247  * - `System::variables_tag`
248  * - Adds: nothing
249  * - Removes: nothing
250  * - Modifies:
251  * - `subcell::Tags::NeighborDataForReconstructionAndRdmpTci<Dim>`
252  * - `evolution::dg::Tags::MortarData`
253  * - `evolution::dg::Tags::MortarNextTemporalId`
254  */
255 template <size_t Dim>
257  template <typename DbTags, typename... InboxTags, typename ArrayIndex,
258  typename ActionList, typename ParallelComponent,
259  typename Metavariables>
261  db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
262  const Parallel::GlobalCache<Metavariables>& /*cache*/,
263  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
264  const ParallelComponent* const /*meta*/) noexcept {
265  static_assert(
266  not Metavariables::local_time_stepping,
267  "DG-subcell does not yet support local time stepping. The "
268  "reconstruction data must be sent using dense output sometimes, and "
269  "not at all other times. However, the data for the RDMP TCI should be "
270  "sent along with the data for reconstruction each time.");
271  const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
272  const auto number_of_expected_messages = element.neighbors().size();
273  if (UNLIKELY(number_of_expected_messages == 0)) {
274  // We have no neighbors, so just continue without doing any work
275  return {std::move(box), Parallel::AlgorithmExecution::Continue};
276  }
277 
278  using ::operator<<;
280  const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
282  FixedHashMap<
283  maximum_number_of_neighbors(Dim), Key,
286  boost::hash<Key>>>& inbox =
288  Metavariables::volume_dim>>(inboxes);
289  const auto& received = inbox.find(current_time_step_id);
290  // Check we have at least some data from correct time, and then check that
291  // we have received all data
292  if (received == inbox.end() or
293  received->second.size() != number_of_expected_messages) {
294  return {std::move(box), Parallel::AlgorithmExecution::Retry};
295  }
296 
297  // Now that we have received all the data, copy it over as needed.
301  boost::hash<Key>>
302  received_data = std::move(inbox[current_time_step_id]);
303  inbox.erase(current_time_step_id);
304  db::mutate<Tags::NeighborDataForReconstructionAndRdmpTci<Dim>,
307  make_not_null(&box),
308  [&current_time_step_id, &element, &received_data](
312  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>*>
313  neighbor_data_ptr,
315  Key, evolution::dg::MortarData<Dim>, boost::hash<Key>>*>
316  mortar_data,
317  const gsl::not_null<
319  mortar_next_time_step_id) noexcept {
320  // Get the next time step id, and also the fluxes data if the neighbor
321  // is doing DG.
322  for (auto& received_mortar_data : received_data) {
323  const auto& mortar_id = received_mortar_data.first;
324  try {
325  mortar_next_time_step_id->at(mortar_id) =
326  std::get<3>(received_mortar_data.second);
327  } catch (std::exception& e) {
328  ERROR("Failed retrieving the MortarId: ("
329  << mortar_id.first << ',' << mortar_id.second
330  << ") from the mortar_next_time_step_id. Got exception: "
331  << e.what());
332  }
333  if (std::get<2>(received_mortar_data.second).has_value()) {
334  mortar_data->at(mortar_id).insert_neighbor_mortar_data(
335  current_time_step_id,
336  std::get<0>(received_mortar_data.second),
337  std::move(*std::get<2>(received_mortar_data.second)));
338  }
339  }
340 
341  ASSERT(neighbor_data_ptr->size() == 1,
342  "Should only have one element in the neighbor data when "
343  "receiving neighbor data");
344  ASSERT(
345  neighbor_data_ptr->count(
346  std::pair{Direction<Dim>::lower_xi(),
347  element.id().external_boundary_id()}) == 1,
348  "The self data for the RDMP TCI should have been inserted but it "
349  "wasn't, despite there being one entry in the neighbor data.");
350  NeighborData& self_neighbor_data = (*neighbor_data_ptr)[std::pair{
353  const size_t number_of_evolved_vars =
354  self_neighbor_data.max_variables_values.size();
355  ASSERT(self_neighbor_data.min_variables_values.size() ==
356  number_of_evolved_vars,
357  "The number of evolved variables for which we have a maximum "
358  "and minimum should be the same, but we have "
359  << number_of_evolved_vars << " for the max and "
360  << self_neighbor_data.min_variables_values.size()
361  << " for the min.");
362 
363  for (const auto& [direction, neighbors_in_direction] :
364  element.neighbors()) {
365  for (const auto& neighbor : neighbors_in_direction) {
366  std::pair directional_element_id{direction, neighbor};
367  ASSERT(neighbor_data_ptr->count(directional_element_id) == 0,
368  "Found neighbor already inserted in direction "
369  << direction << " with ElementId " << neighbor);
370  ASSERT(std::get<1>(received_data[directional_element_id])
371  .has_value(),
372  "Received subcell data message that does not contain any "
373  "actual subcell data for reconstruction.");
374  // Collect the max/min of u(t^n) for the RDMP as we receive data.
375  // This reduces the memory footprint.
376  std::vector<double>& received_neighbor_subcell_data =
377  *std::get<1>(received_data[directional_element_id]);
378  const size_t max_offset = received_neighbor_subcell_data.size() -
379  2 * number_of_evolved_vars;
380  const size_t min_offset = received_neighbor_subcell_data.size() -
381  number_of_evolved_vars;
382  for (size_t var_index = 0; var_index < number_of_evolved_vars;
383  ++var_index) {
384  self_neighbor_data.max_variables_values[var_index] = std::max(
385  self_neighbor_data.max_variables_values[var_index],
386  received_neighbor_subcell_data[max_offset + var_index]);
387  self_neighbor_data.min_variables_values[var_index] = std::min(
388  self_neighbor_data.min_variables_values[var_index],
389  received_neighbor_subcell_data[min_offset + var_index]);
390  }
391  // Copy over the ghost cell data for subcell reconstruction.
392  [[maybe_unused]] const auto insert_result =
393  neighbor_data_ptr->insert(std::pair{
394  directional_element_id,
396  received_neighbor_subcell_data.begin(),
397  std::prev(
398  received_neighbor_subcell_data.end(),
399  2 * static_cast<typename std::iterator_traits<
401  difference_type>(
402  number_of_evolved_vars))}}});
403  ASSERT(insert_result.second,
404  "Failed to insert the neighbor data in direction "
405  << direction << " from neighbor " << neighbor);
406  }
407  }
408  });
409  return {std::move(box), Parallel::AlgorithmExecution::Continue};
410  }
411 };
412 } // namespace evolution::dg::subcell::Actions
evolution::dg::subcell::rdmp_max_min
std::pair< std::vector< double >, std::vector< double > > rdmp_max_min(const Variables< tmpl::list< EvolvedVarsTags... >> &active_grid_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... >> &inactive_grid_evolved_vars, const bool include_inactive_grid) noexcept
get the max and min of each component of the active and inactive variables. If include_inactive_grid ...
Definition: RdmpTci.hpp:125
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
maximum_number_of_neighbors
constexpr size_t maximum_number_of_neighbors(const size_t dim)
Definition: MaxNumberOfNeighbors.hpp:15
ElementId::external_boundary_id
static ElementId< VolumeDim > external_boundary_id() noexcept
Returns an ElementId meant for identifying data on external boundaries, which should never correspond...
evolution::dg::subcell::Actions
Actions for the DG-subcell hybrid solver.
Definition: Actions.hpp:9
orient_variables
Variables< TagsList > orient_variables(const Variables< TagsList > &variables, const Index< VolumeDim > &extents, const OrientationMap< VolumeDim > &orientation_of_neighbor) noexcept
Definition: OrientationMapHelpers.hpp:55
utility
Parallel::AlgorithmExecution::Retry
@ Retry
Temporarily stop executing iterable actions, but try the same action again after receiving data from ...
std::exception
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
std::pair
GlobalCache.hpp
Tags.hpp
vector
MakeArray.hpp
evolution::dg::subcell::Actions::ReceiveDataForReconstruction
Receive the subcell data from our neighbor, and accumulate the data from the relaxed discrete maximum...
Definition: ReconstructionCommunication.hpp:256
iterator
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Index
Definition: Index.hpp:31
Direction::lower_xi
static Direction< VolumeDim > lower_xi() noexcept
Helper functions for creating specific Directions. These are labeled by the logical-coordinate names ...
Spectral.hpp
algorithm
Element
Definition: Element.hpp:29
ElementId< Dim >
ElementId.hpp
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:50
Parallel::AlgorithmExecution::Continue
@ Continue
Leave the algorithm termination flag in its current state.
DataBox.hpp
cstddef
Assert.hpp
std::iterator_traits
array
Index.hpp
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
map
Parallel::AlgorithmExecution
AlgorithmExecution
The possible options for altering the current execution of the algorithm, used in the return type of ...
Definition: AlgorithmMetafunctions.hpp:31
evolution::dg::subcell::Actions::SendDataForReconstruction
Sets the local data from the relaxed discrete maximum principle troubled-cell indicator and sends gho...
Definition: ReconstructionCommunication.hpp:95
TimeStepId
Definition: TimeStepId.hpp:25
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
Parallel::receive_data
void receive_data(Proxy &&proxy, typename ReceiveTag::temporal_id temporal_id, ReceiveDataType &&receive_data, const bool enable_if_disabled=false) noexcept
Send the data args... to the algorithm running on proxy, and tag the message with the identifier temp...
Definition: Invoke.hpp:32
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
evolution::dg::subcell::NeighborData
Holds neighbor data needed for the TCI and reconstruction.
Definition: NeighborData.hpp:24
DirectionMap
Definition: DirectionMap.hpp:15
TimeStepId.hpp
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:382
limits
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:1247
evolution::dg::Tags::MortarData
Data on mortars, indexed by (Direction, ElementId) pairs.
Definition: MortarTags.hpp:27
Tensor.hpp
evolution::dg::subcell::slice_data
DirectionMap< Dim, std::vector< double > > slice_data(const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t number_of_ghost_points, const DirectionMap< Dim, bool > &directions_to_slice) noexcept
Slice the subcell variables needed for neighbors to perform reconstruction.
Definition: SliceData.hpp:45
evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox
The inbox tag for boundary correction communication and DG-subcell ghost zone cells.
Definition: InboxTags.hpp:94
std::optional
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Direction.hpp
evolution::dg::Tags::MortarNextTemporalId
The next temporal id at which to receive data on the specified mortar.
Definition: MortarTags.hpp:78
FixedHashMap
A hash table with a compile-time specified maximum size and ability to efficiently handle perfect has...
Definition: FixedHashMap.hpp:82
Prefixes.hpp
std::unordered_map
TMPL.hpp
Element::neighbors
const Neighbors_t & neighbors() const noexcept
Information about the neighboring Elements.
Definition: Element.hpp:66
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13