InternalMortarDataImpl.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/functional/hash.hpp>
7 #include <cstddef>
8 #include <optional>
9 #include <type_traits>
10 #include <unordered_map>
11 #include <utility>
12 #include <vector>
13 
15 #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 #include "DataStructures/DataVector.hpp"
21 #include "Domain/FaceNormal.hpp"
23 #include "Domain/Structure/DirectionMap.hpp"
26 #include "Domain/Tags.hpp"
27 #include "Domain/TagsTimeDependent.hpp"
28 #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
29 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
30 #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
31 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
32 #include "Evolution/DiscontinuousGalerkin/ProjectToBoundary.hpp"
33 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
35 #include "Time/Tags.hpp"
36 #include "Time/TimeStepId.hpp"
37 #include "Utilities/Gsl.hpp"
38 #include "Utilities/TMPL.hpp"
39 
40 namespace evolution::dg::Actions::detail {
41 template <typename System, size_t Dim, typename BoundaryCorrection,
42  typename TemporaryTags, typename... PackageDataVolumeArgs>
43 void internal_mortar_data_impl(
44  const gsl::not_null<
45  DirectionMap<Dim, std::optional<Variables<tmpl::list<
48  normal_covector_and_magnitude_ptr,
52  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>*>
53  mortar_data_ptr,
54  const BoundaryCorrection& boundary_correction,
55  const Variables<typename System::variables_tag::tags_list>&
56  volume_evolved_vars,
57  const Variables<
58  db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
59  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
60  const Variables<TemporaryTags>& volume_temporaries,
61  const Variables<get_primitive_vars_tags_from_system<System>>* const
62  volume_primitive_variables,
63  const Element<Dim>& element, const Mesh<Dim>& volume_mesh,
64  const std::unordered_map<
66  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>& mortar_meshes,
67  const std::unordered_map<
70  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>& mortar_sizes,
71  const TimeStepId& temporal_id,
73  moving_mesh_map,
74  const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
75  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
76  volume_inverse_jacobian,
77  const PackageDataVolumeArgs&... package_data_volume_args) noexcept {
78  using variables_tags = typename System::variables_tag::tags_list;
79  using flux_variables = typename System::flux_variables;
80  using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
81  tmpl::size_t<Dim>, Frame::Inertial>;
82  using temporary_tags_for_face =
83  typename BoundaryCorrection::dg_package_data_temporary_tags;
84  using primitive_tags_for_face = typename detail::get_primitive_vars<
85  System::has_primitive_and_conservative_vars>::
86  template f<BoundaryCorrection>;
87  using mortar_tags_list = typename BoundaryCorrection::dg_package_field_tags;
88 
89  using dg_package_data_projected_tags =
90  tmpl::append<variables_tags, fluxes_tags, temporary_tags_for_face,
91  primitive_tags_for_face>;
92  Variables<tmpl::remove_duplicates<tmpl::push_back<
93  tmpl::append<dg_package_data_projected_tags,
94  detail::inverse_spatial_metric_tag<System>>,
95  detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>
96  fields_on_face{};
97  std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
98  for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
99  (void)neighbors_in_direction; // unused variable
100  // In order to reduce memory allocations we handle both the upper and
101  // lower neighbors in each direction together since computing the
102  // contributions to the faces is guaranteed to require the same number of
103  // grid points.
104  if (direction.side() == Side::Upper and
105  element.neighbors().count(direction.opposite()) != 0) {
106  continue;
107  }
108 
109  const auto internal_mortars = [&](const Mesh<Dim - 1>& face_mesh,
110  const Direction<Dim>&
111  local_direction) noexcept {
112  // We may not need to bring the volume fluxes or temporaries to the
113  // boundary since that depends on the specific boundary correction we
114  // are using. Silence compilers warnings about them being unused.
115  (void)volume_fluxes;
116  (void)volume_temporaries;
117 
118  // This helper does the following:
119  //
120  // 1. Use a helper function to get data onto the faces. Done either by
121  // slicing (Gauss-Lobatto points) or interpolation (Gauss points).
122  // This is done using the `project_contiguous_data_to_boundary` and
123  // `project_tensors_to_boundary` functions.
124  //
125  // 2. Invoke the boundary correction to get the packaged data. Note
126  // that this is done on the *face* and NOT the mortar.
127  //
128  // 3. Project the packaged data onto the DG mortars (these might need
129  // re-projection onto subcell mortars later).
130 
131  // Perform step 1
133  volume_evolved_vars, volume_mesh,
134  local_direction);
135  if constexpr (tmpl::size<fluxes_tags>::value != 0) {
137  volume_fluxes, volume_mesh,
138  local_direction);
139  }
140  if constexpr (tmpl::size<tmpl::append<
141  temporary_tags_for_face,
142  detail::inverse_spatial_metric_tag<System>>>::value !=
143  0) {
145  tmpl::append<temporary_tags_for_face,
146  detail::inverse_spatial_metric_tag<System>>>(
147  make_not_null(&fields_on_face), volume_temporaries, volume_mesh,
148  local_direction);
149  }
150  if constexpr (System::has_primitive_and_conservative_vars and
151  tmpl::size<primitive_tags_for_face>::value != 0) {
152  ASSERT(volume_primitive_variables != nullptr,
153  "The volume primitive variables are not set even though the "
154  "system has primitive variables.");
155  project_tensors_to_boundary<primitive_tags_for_face>(
156  make_not_null(&fields_on_face), *volume_primitive_variables,
157  volume_mesh, local_direction);
158  } else {
159  (void)volume_primitive_variables;
160  }
161  if (volume_mesh_velocity.has_value()) {
162  if (not face_mesh_velocity.has_value() or
163  (*face_mesh_velocity)[0].size() !=
164  face_mesh.number_of_grid_points()) {
165  face_mesh_velocity =
166  tnsr::I<DataVector, Dim>{face_mesh.number_of_grid_points()};
167  }
168  project_tensor_to_boundary(make_not_null(&*face_mesh_velocity),
169  *volume_mesh_velocity, volume_mesh,
170  local_direction);
171  }
172 
173  // Normalize the normal vectors. We cache the unit normal covector For
174  // flat geometry and static meshes.
175  const bool mesh_is_moving = not moving_mesh_map.is_identity();
176  if (auto& normal_covector_quantity =
177  normal_covector_and_magnitude_ptr->at(local_direction);
178  detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
179  not normal_covector_quantity.has_value()) {
180  if (not normal_covector_quantity.has_value()) {
181  normal_covector_quantity =
182  Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
184  fields_on_face.number_of_grid_points()};
185  }
186  tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
187 
188  for (size_t inertial_index = 0; inertial_index < Dim;
189  ++inertial_index) {
190  volume_unnormalized_normal_covector.get(inertial_index)
191  .set_data_ref(
192  const_cast<double*>( // NOLINT
193  volume_inverse_jacobian
194  .get(local_direction.dimension(), inertial_index)
195  .data()),
196  volume_mesh.number_of_grid_points());
197  }
200  *normal_covector_quantity)),
201  volume_unnormalized_normal_covector, volume_mesh, local_direction);
202 
203  if (local_direction.side() == Side::Lower) {
204  for (auto& normal_covector_component :
206  *normal_covector_quantity)) {
207  normal_covector_component *= -1.0;
208  }
209  }
210 
211  detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
212  make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
213  *normal_covector_quantity)),
215  *normal_covector_quantity)),
216  make_not_null(&fields_on_face),
218  *normal_covector_quantity));
219  }
220 
221  // Perform step 2
222  ASSERT(normal_covector_and_magnitude_ptr->at(local_direction).has_value(),
223  "The magnitude of the normal vector and the unit normal "
224  "covector have not been computed, even though they should "
225  "have been. Direction: "
226  << local_direction);
227 
228  Variables<mortar_tags_list> packaged_data{
229  face_mesh.number_of_grid_points()};
230  // The DataBox is passed in for retrieving the `volume_tags`
231  const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
232  make_not_null(&packaged_data), boundary_correction, fields_on_face,
234  *normal_covector_and_magnitude_ptr->at(local_direction)),
235  face_mesh_velocity, dg_package_data_projected_tags{},
236  package_data_volume_args...);
237  (void)max_abs_char_speed_on_face;
238 
239  // Perform step 3
240  const auto& neighbors_in_local_direction =
241  element.neighbors().at(local_direction);
242  for (const auto& neighbor : neighbors_in_local_direction) {
243  const auto mortar_id = std::make_pair(local_direction, neighbor);
244  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
245  const auto& mortar_size = mortar_sizes.at(mortar_id);
246 
247  // Project the data from the face to the mortar.
248  // Where no projection is necessary we `std::move` the data
249  // directly to avoid a copy. We can't move the data or modify it
250  // in-place when projecting, because in that case the face may
251  // touch two mortars so we need to keep the data around.
252  auto boundary_data_on_mortar =
254  // NOLINTNEXTLINE(bugprone-use-after-move)
255  ? ::dg::project_to_mortar(packaged_data, face_mesh, mortar_mesh,
256  mortar_size)
257  : std::move(packaged_data);
258 
259  // Store the boundary data on this side of the mortar in a way
260  // that is agnostic to the type of boundary correction used. This
261  // currently requires an additional allocation that could be
262  // eliminated either by:
263  //
264  // 1. Having non-owning Variables
265  //
266  // 2. Allow stealing the allocation out of a Variables (and
267  // inserting an allocation).
268  std::vector<double> type_erased_boundary_data_on_mortar{
269  boundary_data_on_mortar.data(),
270  boundary_data_on_mortar.data() + boundary_data_on_mortar.size()};
271  mortar_data_ptr->at(mortar_id).insert_local_mortar_data(
272  temporal_id, face_mesh,
273  std::move(type_erased_boundary_data_on_mortar));
274  }
275  };
276 
277  const Mesh<Dim - 1> face_mesh =
278  volume_mesh.slice_away(direction.dimension());
279 
280  if (fields_on_face.number_of_grid_points() !=
281  face_mesh.number_of_grid_points()) {
282  fields_on_face.initialize(face_mesh.number_of_grid_points());
283  }
284  internal_mortars(face_mesh, direction);
285 
286  if (element.neighbors().count(direction.opposite()) != 0) {
287  if (fields_on_face.number_of_grid_points() !=
288  face_mesh.number_of_grid_points()) {
289  fields_on_face.initialize(face_mesh.number_of_grid_points());
290  }
291  internal_mortars(face_mesh, direction.opposite());
292  }
293  }
294 }
295 
296 template <typename System, size_t Dim, typename BoundaryCorrection,
297  typename DbTagsList, typename... PackageDataVolumeTags>
298 void internal_mortar_data(
299  const gsl::not_null<db::DataBox<DbTagsList>*> box,
300  const BoundaryCorrection& boundary_correction,
301  const Variables<typename System::variables_tag::tags_list>
302  evolved_variables,
303  const Variables<
304  db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
305  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
306  const Variables<
307  typename System::compute_volume_time_derivative_terms::temporary_tags>&
308  temporaries,
309  const Variables<get_primitive_vars_tags_from_system<System>>* const
310  primitive_vars,
311  tmpl::list<PackageDataVolumeTags...> /*meta*/) {
312  db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
314  box,
315  [&boundary_correction,
316  &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
317  &logical_to_inertial_inverse_jacobian = db::get<
319  *box),
320  &mesh = db::get<domain::Tags::Mesh<Dim>>(*box),
321  &mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(*box),
322  &mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box),
323  &mortar_sizes = db::get<Tags::MortarSize<Dim>>(*box),
325  Dim, Frame::Grid, Frame::Inertial>>(*box),
326  &primitive_vars, &temporaries,
327  &time_step_id = db::get<::Tags::TimeStepId>(*box),
328  &volume_fluxes](const auto normal_covector_and_magnitude_ptr,
329  const auto mortar_data_ptr,
330  const auto&... package_data_volume_args) noexcept {
331  detail::internal_mortar_data_impl<System>(
332  normal_covector_and_magnitude_ptr, mortar_data_ptr,
333  boundary_correction, evolved_variables, volume_fluxes, temporaries,
334  primitive_vars, element, mesh, mortar_meshes, mortar_sizes,
335  time_step_id, moving_mesh_map, mesh_velocity,
336  logical_to_inertial_inverse_jacobian, package_data_volume_args...);
337  },
338  db::get<PackageDataVolumeTags>(*box)...);
339 }
340 } // namespace evolution::dg::Actions::detail
Mesh::slice_away
Mesh< Dim - 1 > slice_away(size_t d) const noexcept
Returns a Mesh with dimension d removed (zero-indexed).
Definition: Mesh.cpp:19
utility
Frame::Inertial
Definition: IndexType.hpp:44
domain::CoordinateMapBase
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:45
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
dg::project_to_mortar
Variables< Tags > project_to_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:85
Frame::Grid
Definition: IndexType.hpp:43
std::pair
FaceNormal.hpp
Tags.hpp
vector
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:753
evolution::dg::project_tensor_to_boundary
void project_tensor_to_boundary(const gsl::not_null< Tensor< DataVector, Symm, IndexList > * > face_field, const Tensor< DataVector, Symm, IndexList > &volume_field, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction) noexcept
Projects a tensor to the face.
Definition: ProjectToBoundary.hpp:236
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
CoordinateMap.hpp
Direction< Dim >
Element
Definition: Element.hpp:29
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
ElementId< Dim >
ElementId.hpp
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:44
DataBox.hpp
cstddef
domain::Tags::InverseJacobian< Dim, Frame::Logical, Frame::Inertial >
std::array
evolution::dg::Tags::MagnitudeOfNormal
The magnitude of the unnormalized normal covector to the interface.
Definition: NormalVectorTags.hpp:18
TimeStepId
Definition: TimeStepId.hpp:25
domain::CoordinateMaps::Tags::CoordinateMap
Definition: Tags.hpp:30
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
evolution::dg::project_tensors_to_boundary
void project_tensors_to_boundary(const gsl::not_null< Variables< FaceVarsTagsList > * > face_fields, const Variables< VolumeVarsTagsList > &volume_fields, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction) noexcept
Projects a subset of the tensors in the volume_fields onto the face.
Definition: ProjectToBoundary.hpp:153
Variables.hpp
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:48
DirectionMap
Definition: DirectionMap.hpp:15
TimeStepId.hpp
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
evolution::dg::project_contiguous_data_to_boundary
void project_contiguous_data_to_boundary(const gsl::not_null< Variables< FaceVarsTagsList > * > face_fields, const Variables< VolumeVarsTagsList > &volume_fields, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction) noexcept
Projects a Variables of volume data to a contiguous subset of a boundary Variables
Definition: ProjectToBoundary.hpp:53
evolution::dg::Tags::MortarData
Data on mortars, indexed by (Direction, ElementId) pairs.
Definition: MortarTags.hpp:26
Tensor.hpp
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
optional
evolution::dg::Tags::NormalCovector
The normal covector to the interface.
Definition: NormalVectorTags.hpp:24
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
Mesh.hpp
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Definition: MortarHelpers.cpp:19
Prefixes.hpp
unordered_map
type_traits
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
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13