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 "Evolution/DiscontinuousGalerkin/UsingSubcell.hpp"
27 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
30 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
32 #include "NumericalAlgorithms/Spectral/Projection.hpp"
34 #include "Parallel/AlgorithmMetafunctions.hpp"
35 #include "Parallel/GlobalCache.hpp"
36 #include "Time/Tags.hpp"
37 #include "Time/TimeStepId.hpp"
38 #include "Utilities/Algorithm.hpp"
40 #include "Utilities/Gsl.hpp"
41 #include "Utilities/MakeArray.hpp"
42 #include "Utilities/TMPL.hpp"
43 #include "Utilities/TaggedTuple.hpp"
44 
45 /// \cond
46 namespace evolution::dg::subcell {
47 // We use a forward declaration instead of including a header file to avoid
48 // coupling to the DG-subcell libraries for executables that don't use subcell.
49 template <typename Metavariables, typename DbTagsList>
51  const gsl::not_null<db::DataBox<DbTagsList>*> box,
53  const TimeStepId,
55  maximum_number_of_neighbors(Metavariables::volume_dim),
63  received_temporal_id_and_data) noexcept;
64 } // namespace evolution::dg::subcell
65 /// \endcond
66 
68 namespace detail {
69 template <typename... BoundaryCorrectionTags, typename... Tags,
70  typename BoundaryCorrection>
71 void boundary_correction(
72  const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
73  boundary_corrections_on_mortar,
74  const Variables<tmpl::list<Tags...>>& local_boundary_data,
75  const Variables<tmpl::list<Tags...>>& neighbor_boundary_data,
76  const BoundaryCorrection& boundary_correction,
77  const ::dg::Formulation dg_formulation) noexcept {
78  boundary_correction.dg_boundary_terms(
80  &get<BoundaryCorrectionTags>(*boundary_corrections_on_mortar))...,
81  get<Tags>(local_boundary_data)..., get<Tags>(neighbor_boundary_data)...,
82  dg_formulation);
83 }
84 } // namespace detail
85 
86 /*!
87  * \brief Computes the boundary corrections and lifts them to the volume.
88  *
89  * Given the data from both sides of each mortar, computes the boundary
90  * correction on each mortar and then lifts it into the volume.
91  *
92  * Future additions include:
93  * - boundary conditions, both through ghost cells and by changing the time
94  * derivatives.
95  * - support local time stepping (shouldn't be very difficult)
96  *
97  * When using local time stepping the neighbor sends data at the neighbor's
98  * current temporal id. Along with the boundary data, the next temporal id at
99  * which the neighbor will send data is also sent. This is equal to the
100  * neighbor's `::Tags::Next<::Tags::TimeStepId>`. When inserting into the mortar
101  * data history, we insert the received temporal id, that is, the current time
102  * of the neighbor, along with the boundary correction data.
103  */
104 template <typename Metavariables>
106  using inbox_tags =
108  Metavariables::volume_dim>>;
109  using const_global_cache_tags = tmpl::list<
112 
113  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
114  typename ActionList, typename ParallelComponent>
116  apply(db::DataBox<DbTagsList>& box,
118  const Parallel::GlobalCache<Metavariables>& /*cache*/,
119  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
120  const ParallelComponent* const /*meta*/) noexcept;
121 
122  private:
123  template <typename DbTagsList, typename... InboxTags>
124  static void complete_time_step(
125  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
126 
127  template <typename DbTagsList, typename... InboxTags>
128  static void receive_global_time_stepping(
129  gsl::not_null<db::DataBox<DbTagsList>*> box,
131 
132  template <typename DbTagsList, typename... InboxTags>
133  static void receive_local_time_stepping(
134  gsl::not_null<db::DataBox<DbTagsList>*> box,
136 };
137 
138 template <typename Metavariables>
139 template <typename DbTagsList, typename... InboxTags>
141  const gsl::not_null<db::DataBox<DbTagsList>*> box,
142  const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) noexcept {
143  constexpr size_t volume_dim = Metavariables::system::volume_dim;
144 
145  // Move inbox contents into the DataBox
147  std::map<
148  TimeStepId,
149  FixedHashMap<
150  maximum_number_of_neighbors(volume_dim), Key,
153  boost::hash<Key>>>& inbox =
155  volume_dim>>(*inboxes);
156  const auto& temporal_id = get<::Tags::TimeStepId>(*box);
157  const auto received_temporal_id_and_data = inbox.find(temporal_id);
158  if constexpr (using_subcell_v<Metavariables>) {
159  evolution::dg::subcell::neighbor_reconstructed_face_solution<Metavariables>(
160  box, make_not_null(&*received_temporal_id_and_data));
161  }
162 
163  db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
165  box,
166  [&received_temporal_id_and_data](
168  Key, evolution::dg::MortarData<volume_dim>, boost::hash<Key>>*>
169  mortar_data,
170  const gsl::not_null<
172  mortar_next_time_step_id) noexcept {
173  for (auto& received_mortar_data :
174  received_temporal_id_and_data->second) {
175  const auto& mortar_id = received_mortar_data.first;
176  ASSERT(Metavariables::local_time_stepping
177  ? mortar_next_time_step_id->at(mortar_id) ==
178  received_temporal_id_and_data->first
179  : received_temporal_id_and_data->first ==
180  mortar_data->at(mortar_id).time_step_id(),
181  "Expected to receive mortar data on mortar "
182  << mortar_id << " at time "
183  << mortar_next_time_step_id->at(mortar_id)
184  << " but actually received at time "
185  << received_temporal_id_and_data->first);
186  mortar_next_time_step_id->at(mortar_id) =
187  std::get<3>(received_mortar_data.second);
188  ASSERT(using_subcell_v<Metavariables> or
189  std::get<2>(received_mortar_data.second).has_value(),
190  "Must receive number boundary correction data when not using "
191  "DG-subcell.");
192  if (std::get<2>(received_mortar_data.second).has_value()) {
193  mortar_data->at(mortar_id).insert_neighbor_mortar_data(
194  received_temporal_id_and_data->first,
195  std::get<0>(received_mortar_data.second),
196  std::move(*std::get<2>(received_mortar_data.second)));
197  }
198  }
199  });
200  inbox.erase(received_temporal_id_and_data);
201 }
202 
203 template <typename Metavariables>
204 template <typename DbTagsList, typename... InboxTags>
205 void ApplyBoundaryCorrections<Metavariables>::receive_local_time_stepping(
206  const gsl::not_null<db::DataBox<DbTagsList>*> box,
207  const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) noexcept {
208  constexpr size_t volume_dim = Metavariables::system::volume_dim;
209 
210  using variables_tag = typename Metavariables::system::variables_tag;
211  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
212 
213  // The boundary history coupling computation (which computes the _lifted_
214  // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
215  // using the `NormalDotNumericalFlux` prefix tag. This is because the
216  // returned quantity is more a `dt` quantity than a
217  // `NormalDotNormalDotFlux` since it's been lifted to the volume.
219  std::map<
220  TimeStepId,
221  FixedHashMap<
222  maximum_number_of_neighbors(volume_dim), Key,
225  boost::hash<Key>>>& inbox =
227  volume_dim>>(*inboxes);
228  const auto& local_next_temporal_id =
229  db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
230 
232  volume_dim, typename dt_variables_tag::type>,
234  box, [&inbox, &local_next_temporal_id](
235  const gsl::not_null<
236  std::unordered_map<Key,
240  typename dt_variables_tag::type>,
241  boost::hash<Key>>*>
242  boundary_data_history,
243  const gsl::not_null<
245  mortar_next_time_step_id) noexcept {
246  // Move received boundary data into boundary history.
247  for (auto received_data = inbox.begin();
248  received_data != inbox.end() and
249  received_data->first < local_next_temporal_id;
250  received_data = inbox.erase(received_data)) {
251  const auto& receive_temporal_id = received_data->first;
252  // Loop over all mortars for which we received data at this time
253  for (auto& received_mortar_data : received_data->second) {
254  const auto& mortar_id = received_mortar_data.first;
255  MortarData<Metavariables::volume_dim> neighbor_mortar_data{};
256  // Insert:
257  // - the current TimeStepId of the neighbor
258  // - the current face mesh of the neighbor
259  // - the current boundary correction data of the neighbor
260  ASSERT(std::get<2>(received_mortar_data.second).has_value(),
261  "Did not receive boundary correction data from the "
262  "neighbor\nMortarId: "
263  << mortar_id << "\nTimeStepId: " << receive_temporal_id);
264  ASSERT(
265  mortar_next_time_step_id->at(mortar_id) == receive_temporal_id,
266  "Expected to receive mortar data on mortar "
267  << mortar_id << " at time "
268  << mortar_next_time_step_id->at(mortar_id)
269  << " but actually received at time "
270  << receive_temporal_id);
271  mortar_next_time_step_id->at(mortar_id) =
272  std::get<3>(received_mortar_data.second);
273  neighbor_mortar_data.insert_neighbor_mortar_data(
274  receive_temporal_id, std::get<0>(received_mortar_data.second),
275  std::move(*std::get<2>(received_mortar_data.second)));
276  boundary_data_history->at(mortar_id).remote_insert(
277  receive_temporal_id, std::move(neighbor_mortar_data));
278  }
279  }
280  });
281 }
282 
283 template <typename Metavariables>
284 template <typename DbTagsList, typename... InboxTags>
285 void ApplyBoundaryCorrections<Metavariables>::complete_time_step(
286  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
287  constexpr size_t volume_dim = Metavariables::system::volume_dim;
288 
289  using variables_tag = typename Metavariables::system::variables_tag;
290  using variables_tags = typename variables_tag::tags_list;
291  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
292 
293  // Set up helper lambda that will compute and lift the boundary corrections
294  const Mesh<volume_dim>& volume_mesh =
295  db::get<domain::Tags::Mesh<volume_dim>>(*box);
296  ASSERT(
297  volume_mesh.quadrature() ==
298  make_array<volume_dim>(volume_mesh.quadrature(0)),
299  "Must have isotropic quadrature, but got volume mesh: " << volume_mesh);
300  const bool using_gauss_lobatto_points =
301  volume_mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto;
302 
303  const bool local_time_stepping = Metavariables::local_time_stepping;
304  Scalar<DataVector> volume_det_inv_jacobian{};
305  Scalar<DataVector> volume_det_jacobian{};
306  if (not local_time_stepping and not using_gauss_lobatto_points) {
307  get(volume_det_inv_jacobian)
308  .set_data_ref(make_not_null(&const_cast<DataVector&>(
309  get(db::get<
311  *box)))));
312  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
313  }
314 
315  const auto& mortar_meshes = db::get<Tags::MortarMesh<volume_dim>>(*box);
316  const auto& mortar_sizes = db::get<Tags::MortarSize<volume_dim>>(*box);
317  const ::dg::Formulation dg_formulation =
318  db::get<::dg::Tags::Formulation>(*box);
319 
320  const DirectionMap<volume_dim,
321  std::optional<Variables<tmpl::list<
324  face_normal_covector_and_magnitude =
325  db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>>(
326  *box);
327 
328  const TimeDelta& time_step = db::get<::Tags::TimeStep>(*box);
329  const auto& time_stepper =
330  db::get<typename Metavariables::time_stepper_tag>(*box);
331 
332  const auto compute_and_lift_boundary_corrections =
333  [&dg_formulation, &face_normal_covector_and_magnitude,
334  local_time_stepping, &mortar_meshes, &mortar_sizes, &time_step,
335  &time_stepper, using_gauss_lobatto_points, &volume_det_jacobian,
336  &volume_det_inv_jacobian, &volume_mesh](
337  const auto dt_variables_ptr, const auto variables_ptr,
338  const auto mortar_data_ptr, const auto mortar_data_history_ptr,
339  const auto& boundary_correction) noexcept {
340  using mortar_tags_list = typename std::decay_t<decltype(
341  boundary_correction)>::dg_package_field_tags;
342 
343  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
344  dt_boundary_correction_on_mortar{};
345  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
346  dt_boundary_correction_projected_onto_face{};
347 
349  ElementId<Metavariables::volume_dim>>* mortar_id_ptr =
350  nullptr;
351 
352  const auto compute_correction_coupling =
353  [&boundary_correction, dg_formulation,
354  &dt_boundary_correction_on_mortar,
355  &dt_boundary_correction_projected_onto_face, &dt_variables_ptr,
356  &face_normal_covector_and_magnitude, local_time_stepping,
357  &mortar_id_ptr, &mortar_meshes, &mortar_sizes,
358  using_gauss_lobatto_points, &volume_det_jacobian,
359  &volume_det_inv_jacobian, &volume_mesh](
360  const MortarData<volume_dim>& local_mortar_data,
361  const MortarData<volume_dim>& neighbor_mortar_data) noexcept
363  // Clang thinks we don't need to capture local_time_stepping.
364  (void)local_time_stepping;
365  ASSERT(mortar_id_ptr != nullptr,
366  "Mortar ID pointer should never be nullptr, it must not have "
367  "been set before invoking the lambda.");
368  if (local_time_stepping and not using_gauss_lobatto_points) {
369  // This needs to be updated every call because the Jacobian may be
370  // time-dependent. In the case of time-independent maps and local
371  // time stepping we could first perform the integral on the
372  // boundaries, and then lift to the volume. This is left as a future
373  // optimization.
374  local_mortar_data.get_local_volume_det_inv_jacobian(
375  make_not_null(&volume_det_inv_jacobian));
376  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
377  }
378  const auto& mortar_id = *mortar_id_ptr;
379  const auto& direction = mortar_id.first;
380  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
381 
382  // Extract local and neighbor data, copy into Variables because
383  // we store them in a std::vector for type erasure.
384  const std::pair<Mesh<volume_dim - 1>, std::vector<double>>&
385  local_mesh_and_data = *local_mortar_data.local_mortar_data();
386  const std::pair<Mesh<volume_dim - 1>, std::vector<double>>&
387  neighbor_mesh_and_data =
388  *neighbor_mortar_data.neighbor_mortar_data();
389  Variables<mortar_tags_list> local_data_on_mortar{
390  mortar_mesh.number_of_grid_points()};
391  Variables<mortar_tags_list> neighbor_data_on_mortar{
392  mortar_mesh.number_of_grid_points()};
393  std::copy(std::get<1>(local_mesh_and_data).begin(),
394  std::get<1>(local_mesh_and_data).end(),
395  local_data_on_mortar.data());
396  std::copy(std::get<1>(neighbor_mesh_and_data).begin(),
397  std::get<1>(neighbor_mesh_and_data).end(),
398  neighbor_data_on_mortar.data());
399 
400  // The boundary computations and lifting can be further
401  // optimized by in the h-refinement case having only one
402  // allocation for the face and having the projection from the
403  // mortar to the face be done in place. E.g.
404  // local_data_on_mortar and neighbor_data_on_mortar could be
405  // allocated fewer times, as well as `needs_projection` section
406  // below could do an in-place projection.
407  if (dt_boundary_correction_on_mortar.number_of_grid_points() !=
408  mortar_mesh.number_of_grid_points()) {
409  dt_boundary_correction_on_mortar.initialize(
410  mortar_mesh.number_of_grid_points());
411  }
412 
413  detail::boundary_correction(
414  make_not_null(&dt_boundary_correction_on_mortar),
415  local_data_on_mortar, neighbor_data_on_mortar,
416  boundary_correction, dg_formulation);
417 
418  const std::array<Spectral::MortarSize, volume_dim - 1>& mortar_size =
419  mortar_sizes.at(mortar_id);
420  const Mesh<volume_dim - 1> face_mesh =
421  volume_mesh.slice_away(direction.dimension());
422 
423  auto& dt_boundary_correction =
424  [&dt_boundary_correction_on_mortar,
425  &dt_boundary_correction_projected_onto_face, &face_mesh,
426  &mortar_mesh, &mortar_size]() noexcept
429  mortar_size)) {
430  dt_boundary_correction_projected_onto_face =
431  ::dg::project_from_mortar(dt_boundary_correction_on_mortar,
432  face_mesh, mortar_mesh,
433  mortar_size);
434  return dt_boundary_correction_projected_onto_face;
435  }
436  return dt_boundary_correction_on_mortar;
437  }();
438 
439  Scalar<DataVector> magnitude_of_face_normal{};
440  if (local_time_stepping) {
441  local_mortar_data.get_local_face_normal_magnitude(
442  &magnitude_of_face_normal);
443  } else {
444  ASSERT(face_normal_covector_and_magnitude.count(direction) == 1 and
445  face_normal_covector_and_magnitude.at(direction)
446  .has_value(),
447  "Face normal covector and magnitude not set in "
448  "direction: "
449  << direction);
450  get(magnitude_of_face_normal)
451  .set_data_ref(make_not_null(&const_cast<DataVector&>(
452  get(get<evolution::dg::Tags::MagnitudeOfNormal>(
453  *face_normal_covector_and_magnitude.at(direction))))));
454  }
455 
456  if (using_gauss_lobatto_points) {
457  // The lift_flux function lifts only on the slice, it does not add
458  // the contribution to the volume.
459  ::dg::lift_flux(make_not_null(&dt_boundary_correction),
460  volume_mesh.extents(direction.dimension()),
461  magnitude_of_face_normal);
462  if (local_time_stepping) {
463  return dt_boundary_correction;
464  } else {
465  // Add the flux contribution to the volume data
467  dt_variables_ptr, dt_boundary_correction,
468  volume_mesh.extents(), direction.dimension(),
469  index_to_slice_at(volume_mesh.extents(), direction));
470  }
471  } else {
472  // We are using Gauss points.
473  //
474  // Notes:
475  // - We should really lift both sides simultaneously since this
476  // reduces memory accesses. Lifting all sides at the same time
477  // is unlikely to improve performance since we lift by jumping
478  // through slices. There may also be compatibility issues with
479  // local time stepping.
480  // - If we lift both sides at the same time we first need to deal
481  // with projecting from mortars to the face, then lift off the
482  // faces. With non-owning Variables memory allocations could be
483  // significantly reduced in this code.
484  if (local_time_stepping) {
485  ASSERT(get(volume_det_inv_jacobian).size() > 0,
486  "For local time stepping the volume determinant of the "
487  "inverse Jacobian has not been set.");
488 
489  Scalar<DataVector> face_det_jacobian{};
490  local_mortar_data.get_local_face_det_jacobian(
491  make_not_null(&face_det_jacobian));
492 
493  Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
494  volume_dt_correction{
495  dt_variables_ptr->number_of_grid_points(), 0.0};
497  make_not_null(&volume_dt_correction), volume_det_inv_jacobian,
498  volume_mesh, direction, dt_boundary_correction,
499  magnitude_of_face_normal, face_det_jacobian);
500  return volume_dt_correction;
501  } else {
502  // Project the determinant of the Jacobian to the face. This
503  // could be optimized by caching in the time-independent case.
504  Scalar<DataVector> face_det_jacobian{
505  face_mesh.number_of_grid_points()};
506  const Matrix identity{};
507  auto interpolation_matrices =
508  make_array<volume_dim>(std::cref(identity));
509  const std::pair<Matrix, Matrix>& matrices =
511  volume_mesh.slice_through(direction.dimension()));
512  gsl::at(interpolation_matrices, direction.dimension()) =
513  direction.side() == Side::Upper ? matrices.second
514  : matrices.first;
515  apply_matrices(make_not_null(&get(face_det_jacobian)),
516  interpolation_matrices, get(volume_det_jacobian),
517  volume_mesh.extents());
518 
520  dt_variables_ptr, volume_det_inv_jacobian, volume_mesh,
521  direction, dt_boundary_correction, magnitude_of_face_normal,
522  face_det_jacobian);
523  }
524  }
525 
526  ASSERT(
527  not local_time_stepping,
528  "Must return lifted data when using local time stepping. It's "
529  "likely a missing return statement or incorrect if-else logical "
530  "related to whether or not we are doing local time stepping.");
531  return {};
532  };
533 
534  if constexpr (Metavariables::local_time_stepping) {
535  (void)mortar_data_ptr;
536 
537  for (auto& mortar_id_and_data : *mortar_data_history_ptr) {
538  const auto& mortar_id = mortar_id_and_data.first;
539  auto& mortar_data_history = mortar_id_and_data.second;
540  const auto& direction = mortar_id.first;
541  if (UNLIKELY(mortar_id.second ==
543  ERROR(
544  "Cannot impose boundary conditions on external boundary in "
545  "direction "
546  << direction
547  << " in the ApplyBoundaryCorrections action. Boundary "
548  "conditions are applied in the ComputeTimeDerivative "
549  "action "
550  "instead. You may have unintentionally added external "
551  "mortars in one of the initialization actions.");
552  }
553  mortar_id_ptr = &mortar_id;
554  auto lifted_volume_data = time_stepper.compute_boundary_delta(
555  compute_correction_coupling,
556  make_not_null(&mortar_data_history), time_step);
557 
558  if (using_gauss_lobatto_points) {
559  // Add the flux contribution to the volume data
561  variables_ptr, lifted_volume_data, volume_mesh.extents(),
562  direction.dimension(),
563  index_to_slice_at(volume_mesh.extents(), direction));
564  } else {
565  *variables_ptr += lifted_volume_data;
566  }
567  }
568  } else {
569  (void)time_step;
570  (void)time_stepper;
571  (void)variables_ptr;
572  (void)mortar_data_history_ptr;
573 
574  for (auto& mortar_id_and_data : *mortar_data_ptr) {
575  // Cannot use structured bindings because of a compiler bug where
576  // they cannot be captured by lambdas.
577  const auto& mortar_id = mortar_id_and_data.first;
578  const auto& mortar_data = mortar_id_and_data.second;
579  const auto& direction = mortar_id.first;
580  if (UNLIKELY(mortar_id.second ==
582  ERROR(
583  "Cannot impose boundary conditions on external boundary in "
584  "direction "
585  << direction
586  << " in the ApplyBoundaryCorrections action. Boundary "
587  "conditions are applied in the ComputeTimeDerivative "
588  "action "
589  "instead. You may have unintentionally added external "
590  "mortars in one of the initialization actions.");
591  }
592  mortar_id_ptr = &mortar_id;
593  compute_correction_coupling(mortar_data, mortar_data);
594  // Remove data since it's tagged with the time. In the future we
595  // _might_ be able to reuse allocations, but this optimization
596  // should only be done after profiling.
597  mortar_id_and_data.second.extract();
598  }
599  }
600  };
601 
602  // Now compute the boundary contribution to this element using the helper
603  // lambda
604  const auto& boundary_correction = db::get<
606  *box);
607  using derived_boundary_corrections =
608  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
609 
610  static_assert(
611  tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
612  "All createable classes for boundary corrections must be marked "
613  "final.");
614  tmpl::for_each<derived_boundary_corrections>(
615  [&boundary_correction, &box, &compute_and_lift_boundary_corrections](
616  auto derived_correction_v) noexcept {
617  using DerivedCorrection =
618  tmpl::type_from<decltype(derived_correction_v)>;
619  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
620  // Compute internal boundary quantities on the mortar for sides of
621  // the element that have neighbors, i.e. they are not an external
622  // side.
623  db::mutate<dt_variables_tag, variables_tag,
626  volume_dim, typename dt_variables_tag::type>>(
627  box, compute_and_lift_boundary_corrections,
628  dynamic_cast<const DerivedCorrection&>(boundary_correction));
629  }
630  });
631 }
632 
633 template <typename Metavariables>
634 template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
635  typename ActionList, typename ParallelComponent>
637 ApplyBoundaryCorrections<Metavariables>::apply(
638  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
639  const Parallel::GlobalCache<Metavariables>& /*cache*/,
640  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
641  const ParallelComponent* const /*meta*/) noexcept {
642  constexpr size_t volume_dim = Metavariables::system::volume_dim;
643  const Element<volume_dim>& element =
644  db::get<domain::Tags::Element<volume_dim>>(box);
645 
646  if (UNLIKELY(element.number_of_neighbors() == 0)) {
647  // We have no neighbors, yay!
648  return {std::move(box), Parallel::AlgorithmExecution::Continue};
649  }
650 
651  if constexpr (not Metavariables::local_time_stepping) {
652  const TimeStepId& temporal_id = get<::Tags::TimeStepId>(box);
653  const auto& inbox =
655  Metavariables::volume_dim>>(inboxes);
656  const auto received_temporal_id_and_data = inbox.find(temporal_id);
657  if (received_temporal_id_and_data == inbox.end()) {
658  return {std::move(box), Parallel::AlgorithmExecution::Retry};
659  }
660  const auto& received_neighbor_data = received_temporal_id_and_data->second;
661  for (const auto& [direction, neighbors] : element.neighbors()) {
662  for (const auto& neighbor : neighbors) {
663  const auto neighbor_received =
664  received_neighbor_data.find(std::pair{direction, neighbor});
665  if (neighbor_received == received_neighbor_data.end()) {
666  return {std::move(box), Parallel::AlgorithmExecution::Retry};
667  }
668  }
669  }
670 
671  receive_global_time_stepping(make_not_null(&box), make_not_null(&inboxes));
672  } else {
673  const auto& inbox =
675  Metavariables::volume_dim>>(inboxes);
676 
677  const auto& local_next_temporal_id =
678  db::get<::Tags::Next<::Tags::TimeStepId>>(box);
679  const auto& mortars_next_temporal_id = db::get<
681  box);
682 
683  for (const auto& [mortar_id, mortar_next_temporal_id] :
684  mortars_next_temporal_id) {
685  // If on an external boundary
686  if (UNLIKELY(mortar_id.second ==
688  continue;
689  }
690  auto next_temporal_id = mortar_next_temporal_id;
691  while (next_temporal_id < local_next_temporal_id) {
692  const auto temporal_received = inbox.find(next_temporal_id);
693  if (temporal_received == inbox.end()) {
694  return {std::move(box), Parallel::AlgorithmExecution::Retry};
695  }
696  const auto mortar_received = temporal_received->second.find(mortar_id);
697  if (mortar_received == temporal_received->second.end()) {
698  return {std::move(box), Parallel::AlgorithmExecution::Retry};
699  }
700  next_temporal_id = std::get<3>(mortar_received->second);
701  }
702  }
703 
704  receive_local_time_stepping(make_not_null(&box), make_not_null(&inboxes));
705  }
706 
707  complete_time_step(make_not_null(&box));
708  return {std::move(box), Parallel::AlgorithmExecution::Continue};
709 }
710 } // 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
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....
Parallel::AlgorithmExecution::Retry
@ Retry
Temporarily stop executing iterable actions, but try the same action again after receiving data from ...
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
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
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:791
Spectral.hpp
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
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:51
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:641
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
dg::Tags::Formulation
The DG formulation to use.
Definition: Formulation.hpp:27
Parallel::AlgorithmExecution::Continue
@ Continue
Leave the algorithm termination flag in its current state.
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
Spectral::MortarSize
ChildSize MortarSize
The portion of an element covered by a mortar.
Definition: Projection.hpp:23
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:67
TimeDelta
Definition: Time.hpp:88
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.
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
Parallel::AlgorithmExecution
AlgorithmExecution
The possible options for altering the current execution of the algorithm, used in the return type of ...
Definition: AlgorithmMetafunctions.hpp:31
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:85
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:227
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:49
DirectionMap
Definition: DirectionMap.hpp:15
TimeStepId.hpp
evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.
Definition: Actions.hpp:6
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
evolution::dg::subcell::neighbor_reconstructed_face_solution
void neighbor_reconstructed_face_solution(const gsl::not_null< db::DataBox< DbTagsList > * > box, const gsl::not_null< std::pair< const TimeStepId, FixedHashMap< maximum_number_of_neighbors(Metavariables::volume_dim), std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim >>, std::tuple< Mesh< Metavariables::volume_dim - 1 >, std::optional< std::vector< double >>, std::optional< std::vector< double >>, ::TimeStepId >, boost::hash< std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim >>>>> * > received_temporal_id_and_data) noexcept
Invoked in directions where the neighbor is doing subcell, this function computes the neighbor data o...
Definition: NeighborReconstructedFaceSolution.hpp:65
Gsl.hpp
std::begin
T begin(T... args)
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
evolution::dg::Actions::ApplyBoundaryCorrections
Computes the boundary corrections and lifts them to the volume.
Definition: ApplyBoundaryCorrections.hpp:105
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:196
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::slice_away
Mesh< Dim - 1 > slice_away(size_t d) const noexcept
Returns a Mesh with dimension d removed (zero-indexed).
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:235
Mesh::extents
const Index< Dim > & extents() const noexcept
The number of grid points in each dimension of the grid.
Definition: Mesh.hpp:148
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
Prefixes.hpp
std::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
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