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 "NumericalAlgorithms/Spectral/Projection.hpp"
36 #include "Time/Tags.hpp"
37 #include "Time/TimeStepId.hpp"
38 #include "Utilities/Gsl.hpp"
39 #include "Utilities/TMPL.hpp"
40 
41 namespace evolution::dg::Actions::detail {
42 template <typename System, size_t Dim, typename BoundaryCorrection,
43  typename TemporaryTags, typename... PackageDataVolumeArgs>
44 void internal_mortar_data_impl(
45  const gsl::not_null<
46  DirectionMap<Dim, std::optional<Variables<tmpl::list<
49  normal_covector_and_magnitude_ptr,
53  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>*>
54  mortar_data_ptr,
55  const BoundaryCorrection& boundary_correction,
56  const Variables<typename System::variables_tag::tags_list>&
57  volume_evolved_vars,
58  const Variables<
59  db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
60  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
61  const Variables<TemporaryTags>& volume_temporaries,
62  const Variables<get_primitive_vars_tags_from_system<System>>* const
63  volume_primitive_variables,
64  const Element<Dim>& element, const Mesh<Dim>& volume_mesh,
65  const std::unordered_map<
67  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>& mortar_meshes,
68  const std::unordered_map<
71  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>& mortar_sizes,
72  const TimeStepId& temporal_id,
74  moving_mesh_map,
75  const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
76  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
77  volume_inverse_jacobian,
78  const PackageDataVolumeArgs&... package_data_volume_args) noexcept {
79  using variables_tags = typename System::variables_tag::tags_list;
80  using flux_variables = typename System::flux_variables;
81  using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
82  tmpl::size_t<Dim>, Frame::Inertial>;
83  using temporary_tags_for_face =
84  typename BoundaryCorrection::dg_package_data_temporary_tags;
85  using primitive_tags_for_face = typename detail::get_primitive_vars<
86  System::has_primitive_and_conservative_vars>::
87  template f<BoundaryCorrection>;
88  using mortar_tags_list = typename BoundaryCorrection::dg_package_field_tags;
89 
90  using dg_package_data_projected_tags =
91  tmpl::append<variables_tags, fluxes_tags, temporary_tags_for_face,
92  primitive_tags_for_face>;
93  Variables<tmpl::remove_duplicates<tmpl::push_back<
94  tmpl::append<dg_package_data_projected_tags,
95  detail::inverse_spatial_metric_tag<System>>,
96  detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>
97  fields_on_face{};
98  std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
99  for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
100  (void)neighbors_in_direction; // unused variable
101  // In order to reduce memory allocations we handle both the upper and
102  // lower neighbors in each direction together since computing the
103  // contributions to the faces is guaranteed to require the same number of
104  // grid points.
105  if (direction.side() == Side::Upper and
106  element.neighbors().count(direction.opposite()) != 0) {
107  continue;
108  }
109 
110  const auto internal_mortars = [&](const Mesh<Dim - 1>& face_mesh,
111  const Direction<Dim>&
112  local_direction) noexcept {
113  // We may not need to bring the volume fluxes or temporaries to the
114  // boundary since that depends on the specific boundary correction we
115  // are using. Silence compilers warnings about them being unused.
116  (void)volume_fluxes;
117  (void)volume_temporaries;
118 
119  // This helper does the following:
120  //
121  // 1. Use a helper function to get data onto the faces. Done either by
122  // slicing (Gauss-Lobatto points) or interpolation (Gauss points).
123  // This is done using the `project_contiguous_data_to_boundary` and
124  // `project_tensors_to_boundary` functions.
125  //
126  // 2. Invoke the boundary correction to get the packaged data. Note
127  // that this is done on the *face* and NOT the mortar.
128  //
129  // 3. Project the packaged data onto the DG mortars (these might need
130  // re-projection onto subcell mortars later).
131 
132  // Perform step 1
134  volume_evolved_vars, volume_mesh,
135  local_direction);
136  if constexpr (tmpl::size<fluxes_tags>::value != 0) {
138  volume_fluxes, volume_mesh,
139  local_direction);
140  }
141  if constexpr (tmpl::size<tmpl::append<
142  temporary_tags_for_face,
143  detail::inverse_spatial_metric_tag<System>>>::value !=
144  0) {
146  tmpl::append<temporary_tags_for_face,
147  detail::inverse_spatial_metric_tag<System>>>(
148  make_not_null(&fields_on_face), volume_temporaries, volume_mesh,
149  local_direction);
150  }
151  if constexpr (System::has_primitive_and_conservative_vars and
152  tmpl::size<primitive_tags_for_face>::value != 0) {
153  ASSERT(volume_primitive_variables != nullptr,
154  "The volume primitive variables are not set even though the "
155  "system has primitive variables.");
156  project_tensors_to_boundary<primitive_tags_for_face>(
157  make_not_null(&fields_on_face), *volume_primitive_variables,
158  volume_mesh, local_direction);
159  } else {
160  (void)volume_primitive_variables;
161  }
162  if (volume_mesh_velocity.has_value()) {
163  if (not face_mesh_velocity.has_value() or
164  (*face_mesh_velocity)[0].size() !=
165  face_mesh.number_of_grid_points()) {
166  face_mesh_velocity =
167  tnsr::I<DataVector, Dim>{face_mesh.number_of_grid_points()};
168  }
169  project_tensor_to_boundary(make_not_null(&*face_mesh_velocity),
170  *volume_mesh_velocity, volume_mesh,
171  local_direction);
172  }
173 
174  // Normalize the normal vectors. We cache the unit normal covector For
175  // flat geometry and static meshes.
176  const bool mesh_is_moving = not moving_mesh_map.is_identity();
177  if (auto& normal_covector_quantity =
178  normal_covector_and_magnitude_ptr->at(local_direction);
179  detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
180  not normal_covector_quantity.has_value()) {
181  if (not normal_covector_quantity.has_value()) {
182  normal_covector_quantity =
183  Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
185  fields_on_face.number_of_grid_points()};
186  }
187  tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
188 
189  for (size_t inertial_index = 0; inertial_index < Dim;
190  ++inertial_index) {
191  volume_unnormalized_normal_covector.get(inertial_index)
192  .set_data_ref(
193  const_cast<double*>( // NOLINT
194  volume_inverse_jacobian
195  .get(local_direction.dimension(), inertial_index)
196  .data()),
197  volume_mesh.number_of_grid_points());
198  }
201  *normal_covector_quantity)),
202  volume_unnormalized_normal_covector, volume_mesh, local_direction);
203 
204  if (local_direction.side() == Side::Lower) {
205  for (auto& normal_covector_component :
207  *normal_covector_quantity)) {
208  normal_covector_component *= -1.0;
209  }
210  }
211 
212  detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
213  make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
214  *normal_covector_quantity)),
216  *normal_covector_quantity)),
217  make_not_null(&fields_on_face),
219  *normal_covector_quantity));
220  }
221 
222  // Perform step 2
223  ASSERT(normal_covector_and_magnitude_ptr->at(local_direction).has_value(),
224  "The magnitude of the normal vector and the unit normal "
225  "covector have not been computed, even though they should "
226  "have been. Direction: "
227  << local_direction);
228 
229  Variables<mortar_tags_list> packaged_data{
230  face_mesh.number_of_grid_points()};
231  // The DataBox is passed in for retrieving the `volume_tags`
232  const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
233  make_not_null(&packaged_data), boundary_correction, fields_on_face,
235  *normal_covector_and_magnitude_ptr->at(local_direction)),
236  face_mesh_velocity, dg_package_data_projected_tags{},
237  package_data_volume_args...);
238  (void)max_abs_char_speed_on_face;
239 
240  // Perform step 3
241  const auto& neighbors_in_local_direction =
242  element.neighbors().at(local_direction);
243  for (const auto& neighbor : neighbors_in_local_direction) {
244  const auto mortar_id = std::make_pair(local_direction, neighbor);
245  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
246  const auto& mortar_size = mortar_sizes.at(mortar_id);
247 
248  // Project the data from the face to the mortar.
249  // Where no projection is necessary we `std::move` the data
250  // directly to avoid a copy. We can't move the data or modify it
251  // in-place when projecting, because in that case the face may
252  // touch two mortars so we need to keep the data around.
253  auto boundary_data_on_mortar =
255  // NOLINTNEXTLINE(bugprone-use-after-move)
256  ? ::dg::project_to_mortar(packaged_data, face_mesh, mortar_mesh,
257  mortar_size)
258  : std::move(packaged_data);
259 
260  // Store the boundary data on this side of the mortar in a way
261  // that is agnostic to the type of boundary correction used. This
262  // currently requires an additional allocation that could be
263  // eliminated either by:
264  //
265  // 1. Having non-owning Variables
266  //
267  // 2. Allow stealing the allocation out of a Variables (and
268  // inserting an allocation).
269  std::vector<double> type_erased_boundary_data_on_mortar{
270  boundary_data_on_mortar.data(),
271  boundary_data_on_mortar.data() + boundary_data_on_mortar.size()};
272  mortar_data_ptr->at(mortar_id).insert_local_mortar_data(
273  temporal_id, face_mesh,
274  std::move(type_erased_boundary_data_on_mortar));
275  }
276  };
277 
278  const Mesh<Dim - 1> face_mesh =
279  volume_mesh.slice_away(direction.dimension());
280 
281  if (fields_on_face.number_of_grid_points() !=
282  face_mesh.number_of_grid_points()) {
283  fields_on_face.initialize(face_mesh.number_of_grid_points());
284  }
285  internal_mortars(face_mesh, direction);
286 
287  if (element.neighbors().count(direction.opposite()) != 0) {
288  if (fields_on_face.number_of_grid_points() !=
289  face_mesh.number_of_grid_points()) {
290  fields_on_face.initialize(face_mesh.number_of_grid_points());
291  }
292  internal_mortars(face_mesh, direction.opposite());
293  }
294  }
295 }
296 
297 template <typename System, size_t Dim, typename BoundaryCorrection,
298  typename DbTagsList, typename... PackageDataVolumeTags>
299 void internal_mortar_data(
300  const gsl::not_null<db::DataBox<DbTagsList>*> box,
301  const BoundaryCorrection& boundary_correction,
302  const Variables<typename System::variables_tag::tags_list>
303  evolved_variables,
304  const Variables<
305  db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
306  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
307  const Variables<
308  typename System::compute_volume_time_derivative_terms::temporary_tags>&
309  temporaries,
310  const Variables<get_primitive_vars_tags_from_system<System>>* const
311  primitive_vars,
312  tmpl::list<PackageDataVolumeTags...> /*meta*/) {
313  db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
315  box,
316  [&boundary_correction,
317  &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
318  &logical_to_inertial_inverse_jacobian = db::get<
320  *box),
321  &mesh = db::get<domain::Tags::Mesh<Dim>>(*box),
322  &mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(*box),
323  &mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box),
324  &mortar_sizes = db::get<Tags::MortarSize<Dim>>(*box),
326  Dim, Frame::Grid, Frame::Inertial>>(*box),
327  &primitive_vars, &temporaries,
328  &time_step_id = db::get<::Tags::TimeStepId>(*box),
329  &volume_fluxes](const auto normal_covector_and_magnitude_ptr,
330  const auto mortar_data_ptr,
331  const auto&... package_data_volume_args) noexcept {
332  detail::internal_mortar_data_impl<System>(
333  normal_covector_and_magnitude_ptr, mortar_data_ptr,
334  boundary_correction, evolved_variables, volume_fluxes, temporaries,
335  primitive_vars, element, mesh, mortar_meshes, mortar_sizes,
336  time_step_id, moving_mesh_map, mesh_velocity,
337  logical_to_inertial_inverse_jacobian, package_data_volume_args...);
338  },
339  db::get<PackageDataVolumeTags>(*box)...);
340 }
341 } // namespace evolution::dg::Actions::detail
Spectral::needs_projection
bool needs_projection(const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes) noexcept
Determine whether data needs to be projected between a child mesh and its parent mesh....
utility
Frame::Inertial
Definition: IndexType.hpp:44
domain::CoordinateMapBase
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:45
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:73
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:791
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:233
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:50
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:49
DirectionMap
Definition: DirectionMap.hpp:15
TimeStepId.hpp
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:27
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
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::slice_away
Mesh< Dim - 1 > slice_away(size_t d) const noexcept
Returns a Mesh with dimension d removed (zero-indexed).
Prefixes.hpp
unordered_map
type_traits
dg::mortar_size
MortarSize< Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, size_t dimension, const OrientationMap< Dim > &orientation) noexcept
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13