ApplyBoundaryCorrections.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <map>
8 #include <optional>
9 #include <tuple>
10 #include <type_traits>
11 #include <vector>
12 
15 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
16 #include "Domain/FaceNormal.hpp"
19 #include "Domain/Tags.hpp"
20 #include "Evolution/BoundaryCorrectionTags.hpp"
21 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
22 #include "Evolution/DiscontinuousGalerkin/LiftFromBoundary.hpp"
23 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
24 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
25 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
26 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
32 #include "Parallel/GlobalCache.hpp"
33 #include "Time/Tags.hpp"
34 #include "Time/TimeStepId.hpp"
35 #include "Utilities/Algorithm.hpp"
37 #include "Utilities/Gsl.hpp"
38 #include "Utilities/MakeArray.hpp"
39 #include "Utilities/TMPL.hpp"
40 #include "Utilities/TaggedTuple.hpp"
41 
43 namespace detail {
44 template <typename... BoundaryCorrectionTags, typename... Tags,
45  typename BoundaryCorrection>
46 void boundary_correction(
47  const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
48  boundary_corrections_on_mortar,
49  const Variables<tmpl::list<Tags...>>& local_boundary_data,
50  const Variables<tmpl::list<Tags...>>& neighbor_boundary_data,
51  const BoundaryCorrection& boundary_correction,
52  const ::dg::Formulation dg_formulation) noexcept {
53  boundary_correction.dg_boundary_terms(
55  &get<BoundaryCorrectionTags>(*boundary_corrections_on_mortar))...,
56  get<Tags>(local_boundary_data)..., get<Tags>(neighbor_boundary_data)...,
57  dg_formulation);
58 }
59 } // namespace detail
60 
61 /*!
62  * \brief Computes the boundary corrections and lifts them to the volume.
63  *
64  * Given the data from both sides of each mortar, computes the boundary
65  * correction on each mortar and then lifts it into the volume.
66  *
67  * Future additions include:
68  * - boundary conditions, both through ghost cells and by changing the time
69  * derivatives.
70  * - support local time stepping (shouldn't be very difficult)
71  *
72  * When using local time stepping the neighbor sends data at the neighbor's
73  * current temporal id. Along with the boundary data, the next temporal id at
74  * which the neighbor will send data is also sent. This is equal to the
75  * neighbor's `::Tags::Next<::Tags::TimeStepId>`. When inserting into the mortar
76  * data history, we insert the received temporal id, that is, the current time
77  * of the neighbor, along with the boundary correction data.
78  */
79 template <typename Metavariables>
81  using inbox_tags =
83  Metavariables::volume_dim>>;
84  using const_global_cache_tags = tmpl::list<
87 
88  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
89  typename ActionList, typename ParallelComponent>
90  static std::tuple<db::DataBox<DbTagsList>&&> apply(
91  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
92  const Parallel::GlobalCache<Metavariables>& /*cache*/,
93  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
94  const ParallelComponent* const /*meta*/) noexcept {
95  constexpr size_t volume_dim = Metavariables::system::volume_dim;
96 
98  .number_of_neighbors() == 0)) {
99  // We have no neighbors, yay!
100  return {std::move(box)};
101  }
102 
103  if constexpr (Metavariables::local_time_stepping) {
104  receive_local_time_stepping(make_not_null(&box), make_not_null(&inboxes));
105  } else {
106  receive_global_time_stepping(make_not_null(&box),
107  make_not_null(&inboxes));
108  }
109  complete_time_step(make_not_null(&box));
110  return {std::move(box)};
111  }
112 
113  template <typename DbTags, typename... InboxTags, typename ArrayIndex>
114  static bool is_ready(const db::DataBox<DbTags>& box,
115  const tuples::TaggedTuple<InboxTags...>& inboxes,
116  const Parallel::GlobalCache<Metavariables>& /*cache*/,
117  const ArrayIndex& /*array_index*/) noexcept;
118 
119  private:
120  template <typename DbTagsList, typename... InboxTags>
121  static void complete_time_step(
122  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
123 
124  template <typename DbTagsList, typename... InboxTags>
125  static void receive_global_time_stepping(
126  gsl::not_null<db::DataBox<DbTagsList>*> box,
128 
129  template <typename DbTagsList, typename... InboxTags>
130  static void receive_local_time_stepping(
131  gsl::not_null<db::DataBox<DbTagsList>*> box,
133 };
134 
135 template <typename Metavariables>
136 template <typename DbTagsList, typename... InboxTags>
138  const gsl::not_null<db::DataBox<DbTagsList>*> box,
139  const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) noexcept {
140  constexpr size_t volume_dim = Metavariables::system::volume_dim;
141 
142  // Move inbox contents into the DataBox
144  std::map<
145  TimeStepId,
146  FixedHashMap<
147  maximum_number_of_neighbors(volume_dim), Key,
150  boost::hash<Key>>>& inbox =
152  volume_dim>>(*inboxes);
153  const auto& temporal_id = get<::Tags::TimeStepId>(*box);
154  const auto received_temporal_id_and_data = inbox.find(temporal_id);
155  db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
157  box,
158  [&received_temporal_id_and_data](
160  Key, evolution::dg::MortarData<volume_dim>, boost::hash<Key>>*>
161  mortar_data,
162  const gsl::not_null<
164  mortar_next_time_step_id) noexcept {
165  for (auto& received_mortar_data :
166  received_temporal_id_and_data->second) {
167  const auto& mortar_id = received_mortar_data.first;
168  ASSERT(mortar_next_time_step_id->at(mortar_id) ==
169  received_temporal_id_and_data->first,
170  "Expected to receive mortar data on mortar "
171  << mortar_id << " at time "
172  << mortar_next_time_step_id->at(mortar_id)
173  << " but actually received at time "
174  << received_temporal_id_and_data->first);
175  mortar_next_time_step_id->at(mortar_id) =
176  std::get<3>(received_mortar_data.second);
177  mortar_data->at(mortar_id).insert_neighbor_mortar_data(
178  received_temporal_id_and_data->first,
179  std::get<0>(received_mortar_data.second),
180  std::move(*std::get<2>(received_mortar_data.second)));
181  }
182  });
183  inbox.erase(received_temporal_id_and_data);
184 }
185 
186 template <typename Metavariables>
187 template <typename DbTagsList, typename... InboxTags>
188 void ApplyBoundaryCorrections<Metavariables>::receive_local_time_stepping(
189  const gsl::not_null<db::DataBox<DbTagsList>*> box,
190  const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) noexcept {
191  constexpr size_t volume_dim = Metavariables::system::volume_dim;
192 
193  using variables_tag = typename Metavariables::system::variables_tag;
194  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
195 
196  // The boundary history coupling computation (which computes the _lifted_
197  // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
198  // using the `NormalDotNumericalFlux` prefix tag. This is because the
199  // returned quantity is more a `dt` quantity than a
200  // `NormalDotNormalDotFlux` since it's been lifted to the volume.
202  std::map<
203  TimeStepId,
204  FixedHashMap<
205  maximum_number_of_neighbors(volume_dim), Key,
208  boost::hash<Key>>>& inbox =
210  volume_dim>>(*inboxes);
211  const auto& local_next_temporal_id =
212  db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
213 
215  volume_dim, typename dt_variables_tag::type>,
217  box, [&inbox, &local_next_temporal_id](
218  const gsl::not_null<
219  std::unordered_map<Key,
223  typename dt_variables_tag::type>,
224  boost::hash<Key>>*>
225  boundary_data_history,
226  const gsl::not_null<
228  mortar_next_time_step_id) noexcept {
229  // Move received boundary data into boundary history.
230  for (auto received_data = inbox.begin();
231  received_data != inbox.end() and
232  received_data->first < local_next_temporal_id;
233  received_data = inbox.erase(received_data)) {
234  const auto& receive_temporal_id = received_data->first;
235  // Loop over all mortars for which we received data at this time
236  for (auto& received_mortar_data : received_data->second) {
237  const auto& mortar_id = received_mortar_data.first;
238  MortarData<Metavariables::volume_dim> neighbor_mortar_data{};
239  // Insert:
240  // - the current TimeStepId of the neighbor
241  // - the current face mesh of the neighbor
242  // - the current boundary correction data of the neighbor
243  ASSERT(std::get<2>(received_mortar_data.second).has_value(),
244  "Did not receive boundary correction data from the "
245  "neighbor\nMortarId: "
246  << mortar_id << "\nTimeStepId: " << receive_temporal_id);
247  ASSERT(
248  mortar_next_time_step_id->at(mortar_id) == receive_temporal_id,
249  "Expected to receive mortar data on mortar "
250  << mortar_id << " at time "
251  << mortar_next_time_step_id->at(mortar_id)
252  << " but actually received at time "
253  << receive_temporal_id);
254  mortar_next_time_step_id->at(mortar_id) =
255  std::get<3>(received_mortar_data.second);
256  neighbor_mortar_data.insert_neighbor_mortar_data(
257  receive_temporal_id, std::get<0>(received_mortar_data.second),
258  std::move(*std::get<2>(received_mortar_data.second)));
259  boundary_data_history->at(mortar_id).remote_insert(
260  receive_temporal_id, std::move(neighbor_mortar_data));
261  }
262  }
263  });
264 }
265 
266 template <typename Metavariables>
267 template <typename DbTagsList, typename... InboxTags>
268 void ApplyBoundaryCorrections<Metavariables>::complete_time_step(
269  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
270  constexpr size_t volume_dim = Metavariables::system::volume_dim;
271 
272  using variables_tag = typename Metavariables::system::variables_tag;
273  using variables_tags = typename variables_tag::tags_list;
274  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
275 
276  // Set up helper lambda that will compute and lift the boundary corrections
277  const Mesh<volume_dim>& volume_mesh =
278  db::get<domain::Tags::Mesh<volume_dim>>(*box);
279  ASSERT(
280  volume_mesh.quadrature() ==
281  make_array<volume_dim>(volume_mesh.quadrature(0)),
282  "Must have isotropic quadrature, but got volume mesh: " << volume_mesh);
283  const bool using_gauss_lobatto_points =
284  volume_mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto;
285 
286  const bool local_time_stepping = Metavariables::local_time_stepping;
287  Scalar<DataVector> volume_det_inv_jacobian{};
288  Scalar<DataVector> volume_det_jacobian{};
289  if (not local_time_stepping and not using_gauss_lobatto_points) {
290  get(volume_det_inv_jacobian)
291  .set_data_ref(make_not_null(&const_cast<DataVector&>(
292  get(db::get<
294  *box)))));
295  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
296  }
297 
298  const auto& mortar_meshes = db::get<Tags::MortarMesh<volume_dim>>(*box);
299  const auto& mortar_sizes = db::get<Tags::MortarSize<volume_dim>>(*box);
300  const ::dg::Formulation dg_formulation =
301  db::get<::dg::Tags::Formulation>(*box);
302 
303  const DirectionMap<volume_dim,
304  std::optional<Variables<tmpl::list<
307  face_normal_covector_and_magnitude =
308  db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>>(
309  *box);
310 
311  const TimeDelta& time_step = db::get<::Tags::TimeStep>(*box);
312  const auto& time_stepper =
313  db::get<typename Metavariables::time_stepper_tag>(*box);
314 
315  const auto compute_and_lift_boundary_corrections =
316  [&dg_formulation, &face_normal_covector_and_magnitude,
317  local_time_stepping, &mortar_meshes, &mortar_sizes, &time_step,
318  &time_stepper, using_gauss_lobatto_points, &volume_det_jacobian,
319  &volume_det_inv_jacobian, &volume_mesh](
320  const auto dt_variables_ptr, const auto variables_ptr,
321  const auto mortar_data_ptr, const auto mortar_data_history_ptr,
322  const auto& boundary_correction) noexcept {
323  using mortar_tags_list = typename std::decay_t<decltype(
324  boundary_correction)>::dg_package_field_tags;
325 
326  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
327  dt_boundary_correction_on_mortar{};
328  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
329  dt_boundary_correction_projected_onto_face{};
330 
332  ElementId<Metavariables::volume_dim>>* mortar_id_ptr =
333  nullptr;
334 
335  const auto compute_correction_coupling =
336  [&boundary_correction, dg_formulation,
337  &dt_boundary_correction_on_mortar,
338  &dt_boundary_correction_projected_onto_face, &dt_variables_ptr,
339  &face_normal_covector_and_magnitude, local_time_stepping,
340  &mortar_id_ptr, &mortar_meshes, &mortar_sizes,
341  using_gauss_lobatto_points, &volume_det_jacobian,
342  &volume_det_inv_jacobian, &volume_mesh](
343  const MortarData<volume_dim>& local_mortar_data,
344  const MortarData<volume_dim>& neighbor_mortar_data) noexcept
346  // Clang thinks we don't need to capture local_time_stepping.
347  (void)local_time_stepping;
348  ASSERT(mortar_id_ptr != nullptr,
349  "Mortar ID pointer should never be nullptr, it must not have "
350  "been set before invoking the lambda.");
351  if (local_time_stepping and not using_gauss_lobatto_points) {
352  // This needs to be updated every call because the Jacobian may be
353  // time-dependent. In the case of time-independent maps and local
354  // time stepping we could first perform the integral on the
355  // boundaries, and then lift to the volume. This is left as a future
356  // optimization.
357  local_mortar_data.get_local_volume_det_inv_jacobian(
358  make_not_null(&volume_det_inv_jacobian));
359  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
360  }
361  const auto& mortar_id = *mortar_id_ptr;
362  const auto& direction = mortar_id.first;
363  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
364 
365  // Extract local and neighbor data, copy into Variables because
366  // we store them in a std::vector for type erasure.
367  const std::pair<Mesh<volume_dim - 1>, std::vector<double>>&
368  local_mesh_and_data = *local_mortar_data.local_mortar_data();
369  const std::pair<Mesh<volume_dim - 1>, std::vector<double>>&
370  neighbor_mesh_and_data =
371  *neighbor_mortar_data.neighbor_mortar_data();
372  Variables<mortar_tags_list> local_data_on_mortar{
373  mortar_mesh.number_of_grid_points()};
374  Variables<mortar_tags_list> neighbor_data_on_mortar{
375  mortar_mesh.number_of_grid_points()};
376  std::copy(std::get<1>(local_mesh_and_data).begin(),
377  std::get<1>(local_mesh_and_data).end(),
378  local_data_on_mortar.data());
379  std::copy(std::get<1>(neighbor_mesh_and_data).begin(),
380  std::get<1>(neighbor_mesh_and_data).end(),
381  neighbor_data_on_mortar.data());
382 
383  // The boundary computations and lifting can be further
384  // optimized by in the h-refinement case having only one
385  // allocation for the face and having the projection from the
386  // mortar to the face be done in place. E.g.
387  // local_data_on_mortar and neighbor_data_on_mortar could be
388  // allocated fewer times, as well as `needs_projection` section
389  // below could do an in-place projection.
390  if (dt_boundary_correction_on_mortar.number_of_grid_points() !=
391  mortar_mesh.number_of_grid_points()) {
392  dt_boundary_correction_on_mortar.initialize(
393  mortar_mesh.number_of_grid_points());
394  }
395 
396  detail::boundary_correction(
397  make_not_null(&dt_boundary_correction_on_mortar),
398  local_data_on_mortar, neighbor_data_on_mortar,
399  boundary_correction, dg_formulation);
400 
401  const std::array<Spectral::MortarSize, volume_dim - 1>& mortar_size =
402  mortar_sizes.at(mortar_id);
403  const Mesh<volume_dim - 1> face_mesh =
404  volume_mesh.slice_away(direction.dimension());
405 
406  auto& dt_boundary_correction =
407  [&dt_boundary_correction_on_mortar,
408  &dt_boundary_correction_projected_onto_face, &face_mesh,
409  &mortar_mesh, &mortar_size]() noexcept
411  if (::dg::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
412  dt_boundary_correction_projected_onto_face =
413  ::dg::project_from_mortar(dt_boundary_correction_on_mortar,
414  face_mesh, mortar_mesh,
415  mortar_size);
416  return dt_boundary_correction_projected_onto_face;
417  }
418  return dt_boundary_correction_on_mortar;
419  }();
420 
421  Scalar<DataVector> magnitude_of_face_normal{};
422  if (local_time_stepping) {
423  local_mortar_data.get_local_face_normal_magnitude(
424  &magnitude_of_face_normal);
425  } else {
426  ASSERT(face_normal_covector_and_magnitude.count(direction) == 1 and
427  face_normal_covector_and_magnitude.at(direction)
428  .has_value(),
429  "Face normal covector and magnitude not set in "
430  "direction: "
431  << direction);
432  get(magnitude_of_face_normal)
433  .set_data_ref(make_not_null(&const_cast<DataVector&>(
434  get(get<evolution::dg::Tags::MagnitudeOfNormal>(
435  *face_normal_covector_and_magnitude.at(direction))))));
436  }
437 
438  if (using_gauss_lobatto_points) {
439  // The lift_flux function lifts only on the slice, it does not add
440  // the contribution to the volume.
441  ::dg::lift_flux(make_not_null(&dt_boundary_correction),
442  volume_mesh.extents(direction.dimension()),
443  magnitude_of_face_normal);
444  if (local_time_stepping) {
445  return dt_boundary_correction;
446  } else {
447  // Add the flux contribution to the volume data
449  dt_variables_ptr, dt_boundary_correction,
450  volume_mesh.extents(), direction.dimension(),
451  index_to_slice_at(volume_mesh.extents(), direction));
452  }
453  } else {
454  // We are using Gauss points.
455  //
456  // Notes:
457  // - We should really lift both sides simultaneously since this
458  // reduces memory accesses. Lifting all sides at the same time
459  // is unlikely to improve performance since we lift by jumping
460  // through slices. There may also be compatibility issues with
461  // local time stepping.
462  // - If we lift both sides at the same time we first need to deal
463  // with projecting from mortars to the face, then lift off the
464  // faces. With non-owning Variables memory allocations could be
465  // significantly reduced in this code.
466  if (local_time_stepping) {
467  ASSERT(get(volume_det_inv_jacobian).size() > 0,
468  "For local time stepping the volume determinant of the "
469  "inverse Jacobian has not been set.");
470 
471  Scalar<DataVector> face_det_jacobian{};
472  local_mortar_data.get_local_face_det_jacobian(
473  make_not_null(&face_det_jacobian));
474 
475  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
476  volume_dt_correction{
477  dt_variables_ptr->number_of_grid_points(), 0.0};
479  make_not_null(&volume_dt_correction), volume_det_inv_jacobian,
480  volume_mesh, direction, dt_boundary_correction,
481  magnitude_of_face_normal, face_det_jacobian);
482  return volume_dt_correction;
483  } else {
484  // Project the determinant of the Jacobian to the face. This
485  // could be optimized by caching in the time-independent case.
486  Scalar<DataVector> face_det_jacobian{
487  face_mesh.number_of_grid_points()};
488  const Matrix identity{};
489  auto interpolation_matrices =
490  make_array<volume_dim>(std::cref(identity));
491  const std::pair<Matrix, Matrix>& matrices =
493  volume_mesh.slice_through(direction.dimension()));
494  gsl::at(interpolation_matrices, direction.dimension()) =
495  direction.side() == Side::Upper ? matrices.second
496  : matrices.first;
497  apply_matrices(make_not_null(&get(face_det_jacobian)),
498  interpolation_matrices, get(volume_det_jacobian),
499  volume_mesh.extents());
500 
502  dt_variables_ptr, volume_det_inv_jacobian, volume_mesh,
503  direction, dt_boundary_correction, magnitude_of_face_normal,
504  face_det_jacobian);
505  }
506  }
507 
508  ASSERT(
509  not local_time_stepping,
510  "Must return lifted data when using local time stepping. It's "
511  "likely a missing return statement or incorrect if-else logical "
512  "related to whether or not we are doing local time stepping.");
513  return {};
514  };
515 
516  if constexpr (Metavariables::local_time_stepping) {
517  (void)mortar_data_ptr;
518 
519  for (auto& mortar_id_and_data : *mortar_data_history_ptr) {
520  const auto& mortar_id = mortar_id_and_data.first;
521  auto& mortar_data_history = mortar_id_and_data.second;
522  const auto& direction = mortar_id.first;
523  if (UNLIKELY(mortar_id.second ==
525  ERROR(
526  "Cannot impose boundary conditions on external boundary in "
527  "direction "
528  << direction
529  << " in the ApplyBoundaryCorrections action. Boundary "
530  "conditions are applied in the ComputeTimeDerivative "
531  "action "
532  "instead. You may have unintentionally added external "
533  "mortars in one of the initialization actions.");
534  }
535  mortar_id_ptr = &mortar_id;
536  auto lifted_volume_data = time_stepper.compute_boundary_delta(
537  compute_correction_coupling,
538  make_not_null(&mortar_data_history), time_step);
539 
540  if (using_gauss_lobatto_points) {
541  // Add the flux contribution to the volume data
543  variables_ptr, lifted_volume_data, volume_mesh.extents(),
544  direction.dimension(),
545  index_to_slice_at(volume_mesh.extents(), direction));
546  } else {
547  *variables_ptr += lifted_volume_data;
548  }
549  }
550  } else {
551  (void)time_step;
552  (void)time_stepper;
553  (void)variables_ptr;
554  (void)mortar_data_history_ptr;
555 
556  for (auto& mortar_id_and_data : *mortar_data_ptr) {
557  // Cannot use structured bindings because of a compiler bug where
558  // they cannot be captured by lambdas.
559  const auto& mortar_id = mortar_id_and_data.first;
560  const auto& mortar_data = mortar_id_and_data.second;
561  const auto& direction = mortar_id.first;
562  if (UNLIKELY(mortar_id.second ==
564  ERROR(
565  "Cannot impose boundary conditions on external boundary in "
566  "direction "
567  << direction
568  << " in the ApplyBoundaryCorrections action. Boundary "
569  "conditions are applied in the ComputeTimeDerivative "
570  "action "
571  "instead. You may have unintentionally added external "
572  "mortars in one of the initialization actions.");
573  }
574  mortar_id_ptr = &mortar_id;
575  compute_correction_coupling(mortar_data, mortar_data);
576  }
577  }
578  };
579 
580  // Now compute the boundary contribution to this element using the helper
581  // lambda
582  const auto& boundary_correction = db::get<
584  *box);
585  using derived_boundary_corrections =
586  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
587 
588  static_assert(
589  tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
590  "All createable classes for boundary corrections must be marked "
591  "final.");
592  tmpl::for_each<derived_boundary_corrections>(
593  [&boundary_correction, &box, &compute_and_lift_boundary_corrections](
594  auto derived_correction_v) noexcept {
595  using DerivedCorrection =
596  tmpl::type_from<decltype(derived_correction_v)>;
597  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
598  // Compute internal boundary quantities on the mortar for sides of
599  // the element that have neighbors, i.e. they are not an external
600  // side.
601  db::mutate<dt_variables_tag, variables_tag,
604  volume_dim, typename dt_variables_tag::type>>(
605  box, compute_and_lift_boundary_corrections,
606  dynamic_cast<const DerivedCorrection&>(boundary_correction));
607  }
608  });
609 }
610 
611 template <typename Metavariables>
612 template <typename DbTags, typename... InboxTags, typename ArrayIndex>
613 bool ApplyBoundaryCorrections<Metavariables>::is_ready(
614  const db::DataBox<DbTags>& box,
615  const tuples::TaggedTuple<InboxTags...>& inboxes,
616  const Parallel::GlobalCache<Metavariables>& /*cache*/,
617  const ArrayIndex& /*array_index*/) noexcept {
618  static constexpr size_t volume_dim = Metavariables::volume_dim;
619  const Element<volume_dim>& element =
620  db::get<domain::Tags::Element<volume_dim>>(box);
621 
622  if (UNLIKELY(element.number_of_neighbors() == 0)) {
623  return true;
624  }
625 
626  if constexpr (not Metavariables::local_time_stepping) {
627  const TimeStepId& temporal_id = get<::Tags::TimeStepId>(box);
628  const auto& inbox =
630  Metavariables::volume_dim>>(inboxes);
631  const auto received_temporal_id_and_data = inbox.find(temporal_id);
632  if (received_temporal_id_and_data == inbox.end()) {
633  return false;
634  }
635  const auto& received_neighbor_data = received_temporal_id_and_data->second;
636  for (const auto& [direction, neighbors] : element.neighbors()) {
637  for (const auto& neighbor : neighbors) {
638  const auto neighbor_received =
639  received_neighbor_data.find(std::pair{direction, neighbor});
640  if (neighbor_received == received_neighbor_data.end()) {
641  return false;
642  }
643  }
644  }
645  } else {
646  const auto& inbox =
648  Metavariables::volume_dim>>(inboxes);
649 
650  const auto& local_next_temporal_id =
651  db::get<::Tags::Next<::Tags::TimeStepId>>(box);
652  const auto& mortars_next_temporal_id = db::get<
654  box);
655 
656  for (const auto& [mortar_id, mortar_next_temporal_id] :
657  mortars_next_temporal_id) {
658  // If on an external boundary
659  if (UNLIKELY(mortar_id.second ==
661  continue;
662  }
663  auto next_temporal_id = mortar_next_temporal_id;
664  while (next_temporal_id < local_next_temporal_id) {
665  const auto temporal_received = inbox.find(next_temporal_id);
666  if (temporal_received == inbox.end()) {
667  return false;
668  }
669  const auto mortar_received = temporal_received->second.find(mortar_id);
670  if (mortar_received == temporal_received->second.end()) {
671  return false;
672  }
673  next_temporal_id = std::get<3>(mortar_received->second);
674  }
675  }
676  }
677  return true;
678 }
679 } // namespace evolution::dg::Actions
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
apply_matrices
void apply_matrices(const gsl::not_null< Variables< VariableTags > * > result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:67
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
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
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
index_to_slice_at
size_t index_to_slice_at(const Index< Dim > &extents, const Direction< Dim > &direction, const size_t offset=0) noexcept
Finds the index in the perpendicular dimension of an element boundary.
Definition: IndexToSliceAt.hpp:18
std::pair
GlobalCache.hpp
LiftFlux.hpp
FaceNormal.hpp
Tags.hpp
vector
Error.hpp
MakeArray.hpp
domain::Tags::Element
Definition: Tags.hpp:97
std::size
T size(T... args)
dg::lift_flux
void lift_flux(const gsl::not_null< Variables< tmpl::list< BoundaryCorrectionTags... >> * > boundary_correction_terms, const size_t extent_perpendicular_to_boundary, Scalar< DataVector > magnitude_of_face_normal) noexcept
Lifts the flux contribution from an interface to the volume.
Definition: LiftFlux.hpp:42
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
evolution::dg::lift_boundary_terms_gauss_points
void lift_boundary_terms_gauss_points(const gsl::not_null< Variables< DtTagsList > * > dt_vars, const Scalar< DataVector > &volume_det_inv_jacobian, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction, const Variables< BoundaryCorrectionTagsList > &boundary_corrections, const Scalar< DataVector > &magnitude_of_face_normal, const Scalar< DataVector > &face_det_jacobian) noexcept
Lift the boundary corrections to the volume time derivatives in the specified direction.
Definition: LiftFromBoundary.hpp:72
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:786
Spectral.hpp
Spectral::MortarSize
MortarSize
The portion of an element covered by a mortar.
Definition: Projection.hpp:18
Element
Definition: Element.hpp:29
Spectral::boundary_interpolation_matrices
const std::pair< Matrix, Matrix > & boundary_interpolation_matrices(const Mesh< 1 > &mesh) noexcept
Matrices that interpolate to the lower and upper boundaries of the element.
Definition: Spectral.cpp:396
evolution::Tags::BoundaryCorrection
The boundary correction used for coupling together neighboring cells or applying boundary conditions.
Definition: BoundaryCorrectionTags.hpp:37
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
ElementId.hpp
db::mutate
decltype(auto) mutate(const gsl::not_null< DataBox< TagList > * > box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:636
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:50
dg::Tags::Formulation
The DG formulation to use.
Definition: Formulation.hpp:27
DataBox.hpp
cstddef
add_slice_to_data
void add_slice_to_data(const gsl::not_null< Variables< tmpl::list< VolumeTags... >> * > volume_vars, const Variables< tmpl::list< SliceTags... >> &vars_on_slice, const Index< VolumeDim > &extents, const size_t sliced_dim, const size_t fixed_index) noexcept
Adds data on a codimension 1 slice to a volume quantity. The slice has a constant logical coordinate ...
Definition: SliceVariables.hpp:82
identity
tnsr::Ij< DataType, Dim, Frame::NoFrame > identity(const DataType &used_for_type) noexcept
returns the Identity matrix
Definition: Identity.hpp:14
std::array
evolution::dg::Actions
Actions for using the discontinuous Galerkin to evolve hyperbolic partial differential equations.
Definition: ApplyBoundaryCorrections.hpp:42
TimeDelta
Definition: Time.hpp:88
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
evolution::dg::Tags::MagnitudeOfNormal
The magnitude of the unnormalized normal covector to the interface.
Definition: NormalVectorTags.hpp:18
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
map
TimeStepId
Definition: TimeStepId.hpp:25
std::decay_t
dg::project_from_mortar
Variables< Tags > project_from_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:110
Mesh::slice_through
Mesh< sizeof...(D)> slice_through(D... d) const noexcept
Returns a Mesh with the dimensions d, ... present (zero-indexed).
Definition: Mesh.hpp:218
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
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
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Element::number_of_neighbors
size_t number_of_neighbors() const noexcept
The number of neighbors this element has.
Definition: Element.hpp:69
std::is_final
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
std::begin
T begin(T... args)
evolution::dg::Tags::MortarData
Data on mortars, indexed by (Direction, ElementId) pairs.
Definition: MortarTags.hpp:27
evolution::dg::Actions::ApplyBoundaryCorrections
Computes the boundary corrections and lifts them to the volume.
Definition: ApplyBoundaryCorrections.hpp:80
TimeSteppers::BoundaryHistory
Definition: BoundaryHistory.hpp:31
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
evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox
The inbox tag for boundary correction communication and DG-subcell ghost zone cells.
Definition: InboxTags.hpp:94
optional
Mesh::quadrature
const std::array< Spectral::Quadrature, Dim > & quadrature() const noexcept
The quadrature chosen in each dimension of the grid.
Definition: Mesh.hpp:187
std::end
T end(T... args)
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
evolution::dg::Tags::MortarDataHistory
History of the data on mortars, indexed by (Direction, ElementId) pairs, and used by the linear multi...
Definition: MortarTags.hpp:43
Mesh.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
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:234
Mesh::extents
const Index< Dim > & extents() const noexcept
The number of grid points in each dimension of the grid.
Definition: Mesh.hpp:139
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Definition: MortarHelpers.cpp:19
Prefixes.hpp
std::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
Element::neighbors
const Neighbors_t & neighbors() const noexcept
Information about the neighboring Elements.
Definition: Element.hpp:66
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13