CollectDataForFluxes.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
9 #include "Domain/InterfaceHelpers.hpp"
10 #include "Domain/Tags.hpp"
11 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
12 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
13 #include "NumericalAlgorithms/Spectral/Projection.hpp"
14 #include "Utilities/Gsl.hpp"
15 
16 /// \cond
17 namespace Parallel {
18 template <typename Metavariables>
19 struct GlobalCache;
20 } // namespace Parallel
21 namespace tuples {
22 template <typename... Tags>
23 struct TaggedTuple;
24 } // namespace tuples
25 /// \endcond
26 
27 namespace dg {
28 namespace Actions {
29 
30 /*!
31  * \ingroup ActionsGroup
32  * \ingroup DiscontinuousGalerkinGroup
33  * \brief Collect data that is needed to compute numerical fluxes and store it
34  * on mortars, projecting it if necessary.
35  *
36  * Set the `DirectionsTag` template parameter to
37  * `domain::Tags::InternalDirections` to collect data on internal element
38  * interfaces, so they can be communicated to the element's neighbors in a
39  * subsequent call to `dg::Actions::SendDataForFluxes`.
40  *
41  * Alternatively, set the `DirectionsTag` to
42  * `domain::Tags::BoundaryDirectionsInterior` to collect data on external
43  * element boundaries for imposing boundary conditions through numerical fluxes.
44  * In this case, make sure mortars on external boundaries are initialized e.g.
45  * by `dg::Actions::InitializeMortars`. Also, make sure the field data on
46  * exterior ("ghost") element boundaries (see
47  * `domain::Tags::BoundaryDirectionsExterior`) has been updated to represent the
48  * boundary conditions before invoking this action.
49  *
50  * Design decisions:
51  *
52  * - We assume that all data that is needed on an element to compute numerical
53  * fluxes is also needed on its neighbor. This is reasonable since numerical
54  * fluxes typically satisfy conservation criteria. It is possible that this
55  * assumptions leads to slightly more data being communicated than is necessary,
56  * e.g. in a strong-form scheme with a numerical flux that does not require all
57  * normal-dot-fluxes remotely, but requires them locally to take the difference
58  * to the normal-dot-numerical-fluxes. This overhead is probably negligible,
59  * since the number of projection and communication steps is more relevant than
60  * the amount of data being communicated. This assumption allows us to perform
61  * only one projection to the mortar in each step, instead of projecting local
62  * and remote data separately.
63  * - Terminology: "Boundary data" is data collected on an element interface
64  * that can be projected to a mortar. "Mortar data" is a collection of projected
65  * boundary data from both elements that touch the mortar.
66  */
67 template <typename BoundaryScheme, typename DirectionsTag>
69 
70 /// \cond
71 template <typename BoundaryScheme>
72 struct CollectDataForFluxes<BoundaryScheme, domain::Tags::InternalDirections<
73  BoundaryScheme::volume_dim>> {
74  private:
75  static constexpr size_t volume_dim = BoundaryScheme::volume_dim;
76  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
77  using all_mortar_data_tag =
79  using boundary_data_computer =
80  typename BoundaryScheme::boundary_data_computer;
81 
82  public:
83  template <typename DbTags, typename... InboxTags, typename Metavariables,
84  typename ArrayIndex, typename ActionList,
85  typename ParallelComponent>
86  static std::tuple<db::DataBox<DbTags>&&> apply(
87  db::DataBox<DbTags>& box,
88  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
89  const Parallel::GlobalCache<Metavariables>& /*cache*/,
90  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
91  const ParallelComponent* const /*meta*/) noexcept {
92  // Collect data on element interfaces
93  auto boundary_data_on_interfaces =
94  interface_apply<domain::Tags::InternalDirections<volume_dim>,
95  boundary_data_computer>(box);
96 
97  // Project collected data to all internal mortars and store in DataBox
98  const auto& element = db::get<domain::Tags::Element<volume_dim>>(box);
99  const auto& face_meshes = get<
101  domain::Tags::Mesh<volume_dim - 1>>>(box);
102  const auto& mortar_meshes =
103  get<::Tags::Mortars<domain::Tags::Mesh<volume_dim - 1>, volume_dim>>(
104  box);
105  const auto& mortar_sizes =
106  get<::Tags::Mortars<::Tags::MortarSize<volume_dim - 1>, volume_dim>>(
107  box);
108  const auto& temporal_id = get<temporal_id_tag>(box);
109  for (const auto& direction_and_neighbors : element.neighbors()) {
110  const auto& direction = direction_and_neighbors.first;
111  const auto& face_mesh = face_meshes.at(direction);
112  for (const auto& neighbor : direction_and_neighbors.second) {
113  const auto mortar_id = std::make_pair(direction, neighbor);
114  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
115  const auto& mortar_size = mortar_sizes.at(mortar_id);
116 
117  // Project the data from the face to the mortar.
118  // Where no projection is necessary we `std::move` the data directly to
119  // avoid a copy. We can't move the data or modify it in-place when
120  // projecting, because in that case the face may touch two mortars so we
121  // need to keep the data around.
122  auto boundary_data_on_mortar =
124  ? boundary_data_on_interfaces.at(direction).project_to_mortar(
125  face_mesh, mortar_mesh, mortar_size)
126  : std::move(boundary_data_on_interfaces.at(direction));
127 
128  // Store the boundary data on this side of the mortar
129  db::mutate<all_mortar_data_tag>(
130  make_not_null(&box),
131  [&mortar_id, &temporal_id, &boundary_data_on_mortar ](
133  all_mortar_data) noexcept {
134  all_mortar_data->at(mortar_id).local_insert(
135  temporal_id, std::move(boundary_data_on_mortar));
136  });
137  }
138  }
139  return {std::move(box)};
140  }
141 };
142 
143 template <typename BoundaryScheme>
144 struct CollectDataForFluxes<
145  BoundaryScheme,
146  domain::Tags::BoundaryDirectionsInterior<BoundaryScheme::volume_dim>> {
147  private:
148  static constexpr size_t volume_dim = BoundaryScheme::volume_dim;
149  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
150  using all_mortar_data_tag =
152  using boundary_data_computer =
153  typename BoundaryScheme::boundary_data_computer;
154 
155  public:
156  template <typename DbTags, typename... InboxTags, typename Metavariables,
157  typename ArrayIndex, typename ActionList,
158  typename ParallelComponent>
160  db::DataBox<DbTags>& box,
161  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
162  const Parallel::GlobalCache<Metavariables>& /*cache*/,
163  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
164  const ParallelComponent* const /*meta*/) noexcept {
165  // Collect data on external element boundaries
166  auto interior_boundary_data =
167  interface_apply<domain::Tags::BoundaryDirectionsInterior<volume_dim>,
168  boundary_data_computer>(box);
169  auto exterior_boundary_data =
170  interface_apply<domain::Tags::BoundaryDirectionsExterior<volume_dim>,
171  boundary_data_computer>(box);
172 
173  // Store the boundary data on mortars so the external boundaries are
174  // treated exactly the same as internal boundaries between elements. In
175  // particular, this is important where the
176  // `BoundaryScheme::mortar_data_tag::type` is not just
177  // `SimpleMortarData`, but e.g. keeps track of the mortar history for
178  // local time stepping.
179  const auto& temporal_id = get<temporal_id_tag>(box);
180  for (const auto& direction :
182  const MortarId<volume_dim> mortar_id{
184  db::mutate<all_mortar_data_tag>(
185  make_not_null(&box),
186  [
187  &temporal_id, &mortar_id, &direction, &interior_boundary_data, &
188  exterior_boundary_data
190  all_mortar_data) noexcept {
191  // We don't need to project the boundary data since mortars and
192  // element faces are identical on external boundaries.
193  all_mortar_data->at(mortar_id).local_insert(
194  temporal_id, std::move(interior_boundary_data.at(direction)));
195  all_mortar_data->at(mortar_id).remote_insert(
196  temporal_id, std::move(exterior_boundary_data.at(direction)));
197  });
198  }
199 
200  return {std::move(box)};
201  }
202 };
203 /// \endcond
204 
205 } // namespace Actions
206 } // namespace dg
std::apply
T apply(T... args)
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
Tags::MortarSize
Definition: Tags.hpp:48
Tags.hpp
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
std::tuple
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:296
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ComputeNonconservativeBoundaryFluxes.hpp:23
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...
Definition: ElementId.cpp:73
DataBox.hpp
cstddef
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
dg::Actions::CollectDataForFluxes
Collect data that is needed to compute numerical fluxes and store it on mortars, projecting it if nec...
Definition: CollectDataForFluxes.hpp:68
dg::needs_projection
bool needs_projection(const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:74
Gsl.hpp
domain::Tags::Interface
Tag which is either a SimpleTag for quantities on an interface, base tag to a compute item which acts...
Definition: Tags.hpp:365
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
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Definition: MortarHelpers.cpp:19
dg::mortar_size
std::array< Spectral::MortarSize, Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, const size_t dimension, const OrientationMap< Dim > &orientation) noexcept
Definition: MortarHelpers.cpp:45
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
Tags::Mortars
Definition: Tags.hpp:37