ApplyOperator.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 <string>
9 #include <tuple>
10 #include <type_traits>
11 #include <utility>
12 
14 #include "DataStructures/FixedHashMap.hpp"
15 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
16 #include "Domain/FaceNormal.hpp"
18 #include "Domain/Structure/DirectionMap.hpp"
20 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
21 #include "Domain/Tags.hpp"
22 #include "Domain/Tags/FaceNormal.hpp"
23 #include "Domain/Tags/Faces.hpp"
24 #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
25 #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
26 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
27 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
28 #include "Elliptic/Systems/GetSourcesComputer.hpp"
30 #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
31 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
32 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
34 #include "Parallel/AlgorithmMetafunctions.hpp"
35 #include "Parallel/GlobalCache.hpp"
36 #include "Parallel/InboxInserters.hpp"
37 #include "Parallel/Invoke.hpp"
39 #include "Utilities/GetOutput.hpp"
40 #include "Utilities/Gsl.hpp"
41 #include "Utilities/Literals.hpp"
42 #include "Utilities/TMPL.hpp"
43 #include "Utilities/TaggedTuple.hpp"
44 
45 /// Actions related to elliptic discontinuous Galerkin schemes
47 // The individual actions in this namespace are not exposed publicly because
48 // they don't work on their own. Instead, the public interface (defined below)
49 // exposes them in action lists.
50 namespace detail {
51 
52 // This tag is used to communicate mortar data across neighbors
53 template <size_t Dim, typename TemporalIdTag, typename PrimalFields,
54  typename PrimalFluxes>
55 struct MortarDataInboxTag
57  MortarDataInboxTag<Dim, TemporalIdTag, PrimalFields, PrimalFluxes>> {
58  using temporal_id = typename TemporalIdTag::type;
59  using type = std::map<
60  temporal_id,
63  boost::hash<::dg::MortarId<Dim>>>>;
64 };
65 
66 // Initializes all quantities the DG operator needs on internal and external
67 // faces, as well as the mortars between neighboring elements. Also initializes
68 // the variable-independent background fields in the PDEs.
69 template <typename System, typename BackgroundTag>
70 struct InitializeFacesMortarsAndBackground {
71  private:
72  static constexpr size_t Dim = System::volume_dim;
73  static constexpr bool has_background_fields =
74  not std::is_same_v<typename System::background_fields, tmpl::list<>>;
75  static_assert(
76  not(has_background_fields and std::is_same_v<BackgroundTag, void>),
77  "The system has background fields that need initialization. Supply a "
78  "'BackgroundTag' to 'elliptic::dg::Actions::initialize_operator'.");
79 
82  using InitializeBackground =
84  typename System::background_fields>;
85  using NormalizeFaceNormal =
87 
88  public:
89  using simple_tags = tmpl::append<
90  typename InitializeFacesAndMortars::return_tags,
91  tmpl::conditional_t<has_background_fields,
92  typename InitializeBackground::return_tags,
93  tmpl::list<>>>;
94  using compute_tags = tmpl::list<>;
95 
96  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
97  typename ArrayIndex, typename ActionList,
98  typename ParallelComponent>
99  static std::tuple<db::DataBox<DbTagsList>&&> apply(
100  db::DataBox<DbTagsList>& box,
101  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
102  const Parallel::GlobalCache<Metavariables>& /*cache*/,
103  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
104  const ParallelComponent* const /*meta*/) noexcept {
105  // Initialize faces and mortars
106  db::mutate_apply<typename InitializeFacesAndMortars::return_tags,
107  typename InitializeFacesAndMortars::argument_tags>(
109  db::get<domain::Tags::InitialExtents<Dim>>(box));
110  // Initialize background fields
111  if constexpr (has_background_fields) {
112  db::mutate_apply<typename InitializeBackground::return_tags,
113  typename InitializeBackground::argument_tags>(
115  db::get<BackgroundTag>(box));
116  }
117  // Normalize face normals
118  for (auto& direction : Direction<Dim>::all_directions()) {
119  elliptic::util::mutate_apply_at<
120  domain::make_faces_tags<Dim,
121  typename NormalizeFaceNormal::return_tags>,
122  domain::make_faces_tags<Dim,
123  typename NormalizeFaceNormal::argument_tags>,
124  tmpl::list<>>(NormalizeFaceNormal{}, make_not_null(&box), direction);
125  }
126  return {std::move(box)};
127  }
128 };
129 
130 // Compute auxiliary variables and fluxes from the primal variables, prepare the
131 // local side of all mortars and send the local mortar data to neighbors. Also
132 // handle boundary conditions by preparing the exterior ("ghost") side of
133 // external mortars.
134 template <typename System, bool Linearized, typename TemporalIdTag,
135  typename PrimalFieldsTag, typename PrimalFluxesTag,
136  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
137  typename PrimalMortarFluxesTag,
138  typename FluxesArgsTags =
139  typename System::fluxes_computer::argument_tags,
140  typename SourcesArgsTags = typename elliptic::get_sources_computer<
141  System, Linearized>::argument_tags>
142 struct PrepareAndSendMortarData;
143 
144 template <typename System, bool Linearized, typename TemporalIdTag,
145  typename PrimalFieldsTag, typename PrimalFluxesTag,
146  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
147  typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
148  typename... SourcesArgsTags>
149 struct PrepareAndSendMortarData<
150  System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
151  OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
152  tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
153  private:
154  static constexpr size_t Dim = System::volume_dim;
155  using all_mortar_data_tag = ::Tags::Mortars<
156  elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
157  typename PrimalMortarFieldsTag::tags_list,
158  typename PrimalMortarFluxesTag::tags_list>,
159  Dim>;
160  using mortar_data_inbox_tag =
161  MortarDataInboxTag<Dim, TemporalIdTag,
162  typename PrimalMortarFieldsTag::tags_list,
163  typename PrimalMortarFluxesTag::tags_list>;
164  using BoundaryConditionsBase = typename System::boundary_conditions_base;
165 
166  public:
167  // Request these tags be added to the DataBox by the `SetupDataBox` action. We
168  // don't actually need to initialize them, because the `TemporalIdTag` and the
169  // `PrimalFieldsTag` will be set by other actions before applying the operator
170  // and the remaining tags hold output of the operator.
171  using simple_tags =
172  tmpl::list<TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
173  OperatorAppliedToFieldsTag, all_mortar_data_tag>;
174  using compute_tags = tmpl::list<>;
175 
176  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
177  typename ActionList, typename ParallelComponent>
178  static std::tuple<db::DataBox<DbTagsList>&&> apply(
179  db::DataBox<DbTagsList>& box,
180  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
182  const ElementId<Dim>& element_id, const ActionList /*meta*/,
183  const ParallelComponent* const /*meta*/) noexcept {
184  // Used to retrieve items out of the DataBox to forward to functions
185  const auto get_items = [](const auto&... args) noexcept {
186  return std::forward_as_tuple(args...);
187  };
188  const auto& temporal_id = db::get<TemporalIdTag>(box);
189  const auto& element = db::get<domain::Tags::Element<Dim>>(box);
190  const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
191  const size_t num_points = mesh.number_of_grid_points();
192  const auto& mortar_meshes =
193  db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box);
194  const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
195  const auto& boundary_conditions = domain.blocks()
196  .at(element_id.block_id())
197  .external_boundary_conditions();
198  const auto apply_boundary_condition =
199  [&box, &boundary_conditions, &element_id](
200  const Direction<Dim>& direction,
201  const auto... fields_and_fluxes) noexcept {
202  ASSERT(
203  boundary_conditions.contains(direction),
204  "No boundary condition is available in block "
205  << element_id.block_id() << " in direction " << direction
206  << ". Make sure you are setting up boundary conditions when "
207  "creating the domain.");
208  ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
209  boundary_conditions.at(direction).get()) != nullptr,
210  "The boundary condition in block "
211  << element_id.block_id() << " in direction " << direction
212  << " has an unexpected type. Make sure it derives off the "
213  "'boundary_conditions_base' class set in the system.");
214  const auto& boundary_condition =
215  dynamic_cast<const BoundaryConditionsBase&>(
216  *boundary_conditions.at(direction));
217  elliptic::apply_boundary_condition<Linearized>(
218  boundary_condition, box, direction, fields_and_fluxes...);
219  };
220 
221  // Can't `db::get` the arguments for the boundary conditions within
222  // `db::mutate`, so we extract the data to mutate and move it back in when
223  // we're done.
224  typename PrimalFluxesTag::type primal_fluxes;
225  typename all_mortar_data_tag::type all_mortar_data;
226  db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
227  make_not_null(&box),
228  [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
229  const auto local_all_mortar_data) {
230  primal_fluxes = std::move(*local_primal_fluxes);
231  all_mortar_data = std::move(*local_all_mortar_data);
232  });
233 
234  // Prepare mortar data
235  //
236  // These memory buffers will be discarded when the action returns so we
237  // don't inflate the memory usage of the simulation when the element is
238  // inactive.
239  Variables<typename System::auxiliary_fields> auxiliary_fields_buffer{
240  num_points};
241  Variables<typename System::auxiliary_fluxes> auxiliary_fluxes_buffer{
242  num_points};
243  using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
244  using fluxes_args_volume_tags =
245  typename System::fluxes_computer::volume_tags;
247  fluxes_args_on_faces{};
248  for (const auto& direction : Direction<Dim>::all_directions()) {
249  fluxes_args_on_faces.emplace(
250  direction, elliptic::util::apply_at<
251  domain::make_faces_tags<Dim, fluxes_args_tags,
252  fluxes_args_volume_tags>,
253  fluxes_args_volume_tags>(get_items, box, direction));
254  }
255  elliptic::dg::prepare_mortar_data<System, Linearized>(
256  make_not_null(&auxiliary_fields_buffer),
257  make_not_null(&auxiliary_fluxes_buffer), make_not_null(&primal_fluxes),
258  make_not_null(&all_mortar_data), db::get<PrimalFieldsTag>(box), element,
261  Frame::Inertial>>(box),
265  mortar_meshes,
267  temporal_id, apply_boundary_condition,
268  std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
269  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
270  fluxes_args_on_faces);
271 
272  // Move the mutated data back into the DataBox
273  db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
274  make_not_null(&box),
275  [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
276  const auto local_all_mortar_data) {
277  *local_primal_fluxes = std::move(primal_fluxes);
278  *local_all_mortar_data = std::move(all_mortar_data);
279  });
280 
281  // Send mortar data to neighbors
282  auto& receiver_proxy =
283  Parallel::get_parallel_component<ParallelComponent>(cache);
284  for (const auto& [direction, neighbors] : element.neighbors()) {
285  const size_t dimension = direction.dimension();
286  const auto& orientation = neighbors.orientation();
287  const auto direction_from_neighbor = orientation(direction.opposite());
288  for (const auto& neighbor_id : neighbors) {
289  const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
290  // Make a copy of the local boundary data on the mortar to send to the
291  // neighbor
292  auto remote_boundary_data_on_mortar =
293  get<all_mortar_data_tag>(box).at(mortar_id).local_data(
294  temporal_id);
295  // Reorient the data to the neighbor orientation if necessary
296  if (not orientation.is_aligned()) {
297  remote_boundary_data_on_mortar.orient_on_slice(
298  mortar_meshes.at(mortar_id).extents(), dimension, orientation);
299  }
300  // Send remote data to neighbor
301  Parallel::receive_data<mortar_data_inbox_tag>(
302  receiver_proxy[neighbor_id], temporal_id,
303  std::make_pair(
304  ::dg::MortarId<Dim>{direction_from_neighbor, element.id()},
305  std::move(remote_boundary_data_on_mortar)));
306  } // loop over neighbors in direction
307  } // loop over directions
308 
309  return {std::move(box)};
310  }
311 };
312 
313 // Wait until all mortar data from neighbors is available. Then add boundary
314 // corrections to the primal fluxes, compute their derivatives (i.e. the second
315 // derivatives of the primal variables) and add boundary corrections to the
316 // result.
317 template <typename System, bool Linearized, typename TemporalIdTag,
318  typename PrimalFieldsTag, typename PrimalFluxesTag,
319  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
320  typename PrimalMortarFluxesTag,
321  typename FluxesArgsTags =
322  typename System::fluxes_computer::argument_tags,
323  typename SourcesArgsTags = typename elliptic::get_sources_computer<
324  System, Linearized>::argument_tags>
325 struct ReceiveMortarDataAndApplyOperator;
326 
327 template <typename System, bool Linearized, typename TemporalIdTag,
328  typename PrimalFieldsTag, typename PrimalFluxesTag,
329  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
330  typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
331  typename... SourcesArgsTags>
332 struct ReceiveMortarDataAndApplyOperator<
333  System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
334  OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
335  tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
336  private:
337  static constexpr size_t Dim = System::volume_dim;
338  using all_mortar_data_tag = ::Tags::Mortars<
339  elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
340  typename PrimalMortarFieldsTag::tags_list,
341  typename PrimalMortarFluxesTag::tags_list>,
342  Dim>;
343  using mortar_data_inbox_tag =
344  MortarDataInboxTag<Dim, TemporalIdTag,
345  typename PrimalMortarFieldsTag::tags_list,
346  typename PrimalMortarFluxesTag::tags_list>;
347 
348  public:
349  using const_global_cache_tags =
352  using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
353 
354  template <typename DbTags, typename... InboxTags, typename Metavariables,
355  typename ArrayIndex, typename ActionList,
356  typename ParallelComponent>
358  db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
359  const Parallel::GlobalCache<Metavariables>& /*cache*/,
360  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
361  const ParallelComponent* const /*meta*/) noexcept {
362  const auto& temporal_id = get<TemporalIdTag>(box);
363 
364  if (not ::dg::has_received_from_all_mortars<mortar_data_inbox_tag>(
365  temporal_id, get<domain::Tags::Element<Dim>>(box), inboxes)) {
366  return {std::move(box), Parallel::AlgorithmExecution::Retry};
367  }
368 
369  // Move received "remote" mortar data into the DataBox
370  if (LIKELY(db::get<domain::Tags::Element<Dim>>(box).number_of_neighbors() >
371  0)) {
372  auto received_mortar_data =
373  std::move(tuples::get<mortar_data_inbox_tag>(inboxes)
374  .extract(temporal_id)
375  .mapped());
376  db::mutate<all_mortar_data_tag>(
377  make_not_null(&box), [&received_mortar_data, &temporal_id](
378  const auto all_mortar_data) noexcept {
379  for (auto& [mortar_id, mortar_data] : received_mortar_data) {
380  all_mortar_data->at(mortar_id).remote_insert(
381  temporal_id, std::move(mortar_data));
382  }
383  });
384  }
385 
386  // Apply DG operator
387  db::mutate<OperatorAppliedToFieldsTag, all_mortar_data_tag>(
388  make_not_null(&box),
389  [](const auto&... args) noexcept {
390  elliptic::dg::apply_operator<System, Linearized>(args...);
391  },
392  db::get<PrimalFieldsTag>(box), db::get<PrimalFluxesTag>(box),
393  db::get<domain::Tags::Mesh<Dim>>(box),
395  Frame::Inertial>>(box),
396  db::get<domain::Tags::DetInvJacobian<Frame::Logical, Frame::Inertial>>(
397  box),
400  db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
402  db::get<elliptic::dg::Tags::PenaltyParameter>(box),
403  db::get<elliptic::dg::Tags::Massive>(box), temporal_id,
404  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
405 
406  return {std::move(box), Parallel::AlgorithmExecution::Continue};
407  }
408 };
409 
410 } // namespace detail
411 
412 /*!
413  * \brief Initialize geometric and background quantities for the elliptic DG
414  * operator
415  *
416  * The geometric and background quantities are initialized together because the
417  * geometry depends on the background metric through the normalization of face
418  * normals. Other examples for background fields are curvature quantities
419  * associated with the background metric, or matter sources such as a
420  * mass-density in the XCTS equations. All `System::background_fields` are
421  * retrieved from the `BackgroundTag` together, to enable re-using cached
422  * temporary quantities in the computations. The `variables` function is invoked
423  * on the `BackgroundTag` with the inertial coordinates, the element's `Mesh`
424  * and the element's inverse Jacobian. These arguments allow computing numeric
425  * derivatives, if necessary. The `BackgroundTag` can be set to `void` (default)
426  * if the `System` has no background fields.
427  *
428  * DataBox:
429  * - Uses:
430  * - `domain::Tags::InitialExtents<Dim>`
431  * - `BackgroundTag`
432  * - Adds:
433  * - `::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>`
434  * - `::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>`
435  * - `::Tags::Variables<background_fields>`
436  * - Adds on internal and external faces:
437  * - `domain::Tags::Coordinates<Dim, Frame::Inertial>`
438  * - `::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>`
439  * - `::Tags::Magnitude<domain::Tags::UnnormalizedFaceNormal<Dim>>`
440  * - `::Tags::Variables<background_fields>`
441  *
442  * \note This action relies on the `SetupDataBox` aggregated initialization
443  * mechanism, so `Actions::SetupDataBox` must be present in the `Initialization`
444  * phase action list prior to this action.
445  *
446  * \see elliptic::dg::Actions::apply_operator
447  */
448 template <typename System, typename BackgroundTag = void>
449 using initialize_operator = tmpl::list<
450  detail::InitializeFacesMortarsAndBackground<System, BackgroundTag>>;
451 
452 /*!
453  * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
454  * the `OperatorAppliedToFieldsTag`
455  *
456  * Add this list to the action list of a parallel component to compute the
457  * elliptic DG operator or its linearization. The operator involves a
458  * communication between nearest-neighbor elements. See `elliptic::dg` for
459  * details on the elliptic DG operator. Make sure to add
460  * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
461  * your parallel component so the required DataBox tags are set up before
462  * applying the operator.
463  *
464  * The result of the computation is written to the `OperatorAppliedToFieldsTag`.
465  * Additionally, the primal fluxes are written to the `PrimalFluxesTag` as an
466  * intermediate result. The auxiliary fields and fluxes are discarded to avoid
467  * inflating the memory usage.
468  *
469  * You can specify the `PrimalMortarFieldsTag` and the `PrimalMortarFluxesTag`
470  * to re-use mortar-data memory buffers from other operator applications, for
471  * example when applying the nonlinear and linearized operator. They default to
472  * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
473  * corresponding to these tags are set up in the DataBox.
474  */
475 template <typename System, bool Linearized, typename TemporalIdTag,
476  typename PrimalFieldsTag, typename PrimalFluxesTag,
477  typename OperatorAppliedToFieldsTag,
478  typename PrimalMortarFieldsTag = PrimalFieldsTag,
479  typename PrimalMortarFluxesTag = PrimalFluxesTag>
480 using apply_operator =
481  tmpl::list<detail::PrepareAndSendMortarData<
482  System, Linearized, TemporalIdTag, PrimalFieldsTag,
483  PrimalFluxesTag, OperatorAppliedToFieldsTag,
484  PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
485  detail::ReceiveMortarDataAndApplyOperator<
486  System, Linearized, TemporalIdTag, PrimalFieldsTag,
487  PrimalFluxesTag, OperatorAppliedToFieldsTag,
488  PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
489 
490 /*!
491  * \brief For linear systems, impose inhomogeneous boundary conditions as
492  * contributions to the fixed sources (i.e. the RHS of the equations).
493  *
494  * \see elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
495  */
496 template <
497  typename System, typename FixedSourcesTag,
498  typename FluxesArgsTags = typename System::fluxes_computer::argument_tags,
499  typename SourcesArgsTags = typename System::sources_computer::argument_tags>
501 
502 /// \cond
503 template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
504  typename... SourcesArgsTags>
506  System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
507  tmpl::list<SourcesArgsTags...>> {
508  private:
509  static constexpr size_t Dim = System::volume_dim;
510  using BoundaryConditionsBase = typename System::boundary_conditions_base;
511 
512  public:
513  using const_global_cache_tags =
516 
517  template <typename DbTags, typename... InboxTags, typename Metavariables,
518  typename ActionList, typename ParallelComponent>
519  static std::tuple<db::DataBox<DbTags>&&> apply(
520  db::DataBox<DbTags>& box,
521  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
522  const Parallel::GlobalCache<Metavariables>& /*cache*/,
523  const ElementId<Dim>& element_id, const ActionList /*meta*/,
524  const ParallelComponent* const /*meta*/) noexcept {
525  // Used to retrieve items out of the DataBox to forward to functions
526  const auto get_items = [](const auto&... args) noexcept {
527  return std::forward_as_tuple(args...);
528  };
529  const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
530  const auto& boundary_conditions = domain.blocks()
531  .at(element_id.block_id())
532  .external_boundary_conditions();
533  const auto apply_boundary_condition =
534  [&box, &boundary_conditions, &element_id](
535  const Direction<Dim>& direction,
536  const auto... fields_and_fluxes) noexcept {
537  ASSERT(
538  boundary_conditions.contains(direction),
539  "No boundary condition is available in block "
540  << element_id.block_id() << " in direction " << direction
541  << ". Make sure you are setting up boundary conditions when "
542  "creating the domain.");
543  ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
544  boundary_conditions.at(direction).get()) != nullptr,
545  "The boundary condition in block "
546  << element_id.block_id() << " in direction " << direction
547  << " has an unexpected type. Make sure it derives off the "
548  "'boundary_conditions_base' class set in the system.");
549  const auto& boundary_condition =
550  dynamic_cast<const BoundaryConditionsBase&>(
551  *boundary_conditions.at(direction));
552  elliptic::apply_boundary_condition<false>(
553  boundary_condition, box, direction, fields_and_fluxes...);
554  };
555 
556  // Can't `db::get` the arguments for the boundary conditions within
557  // `db::mutate`, so we extract the data to mutate and move it back in when
558  // we're done.
559  typename FixedSourcesTag::type fixed_sources;
560  db::mutate<FixedSourcesTag>(
561  make_not_null(&box), [&fixed_sources](const auto local_fixed_sources) {
562  fixed_sources = std::move(*local_fixed_sources);
563  });
564 
565  using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
566  using fluxes_args_volume_tags =
567  typename System::fluxes_computer::volume_tags;
569  fluxes_args_on_faces{};
570  for (const auto& direction : Direction<Dim>::all_directions()) {
571  fluxes_args_on_faces.emplace(
572  direction, elliptic::util::apply_at<
573  domain::make_faces_tags<Dim, fluxes_args_tags,
574  fluxes_args_volume_tags>,
575  fluxes_args_volume_tags>(get_items, box, direction));
576  }
577 
578  elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source<System>(
579  make_not_null(&fixed_sources), db::get<domain::Tags::Element<Dim>>(box),
582  Frame::Inertial>>(box),
584  box),
590  db::get<elliptic::dg::Tags::PenaltyParameter>(box),
591  db::get<elliptic::dg::Tags::Massive>(box), apply_boundary_condition,
592  std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
593  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
594  fluxes_args_on_faces);
595 
596  // Move the mutated data back into the DataBox
597  db::mutate<FixedSourcesTag>(
598  make_not_null(&box), [&fixed_sources](const auto local_fixed_sources) {
599  *local_fixed_sources = std::move(fixed_sources);
600  });
601  return {std::move(box)};
602  }
603 };
604 /// \endcond
605 
606 } // namespace elliptic::dg::Actions
elliptic::dg::Actions
Actions related to elliptic discontinuous Galerkin schemes.
Definition: ApplyOperator.hpp:46
maximum_number_of_neighbors
constexpr size_t maximum_number_of_neighbors(const size_t dim)
Definition: MaxNumberOfNeighbors.hpp:15
elliptic::dg::Tags::Massive
Whether or not to multiply the DG operator with the mass matrix. Massive DG operators can be easier t...
Definition: Tags.hpp:66
utility
Frame::Inertial
Definition: IndexType.hpp:44
Parallel::AlgorithmExecution::Retry
@ Retry
Temporarily stop executing iterable actions, but try the same action again after receiving data from ...
elliptic::dg::NormalizeFaceNormal
Normalize face normals, possibly using a background metric. Set InvMetricTag to void to normalize fac...
Definition: Initialization.hpp:143
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
Literals.hpp
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: ApplyOperator.hpp:500
std::pair
GlobalCache.hpp
Tags::MortarSize
Definition: Tags.hpp:48
FaceNormal.hpp
Tags.hpp
domain::Tags::Element
Definition: Tags.hpp:97
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
FixedHashMap::emplace
std::pair< iterator, bool > emplace(Args &&... args) noexcept
Inserts the element if it does not exists.
Definition: FixedHashMap.hpp:141
elliptic::dg::InitializeFacesAndMortars
Initialize the geometry on faces and mortars for the elliptic DG operator.
Definition: Initialization.hpp:74
Direction
Definition: Direction.hpp:23
domain::Tags::Faces
The Tag on element faces.
Definition: Faces.hpp:25
ElementId< Dim >
elliptic::dg::InitializeBackground
Initialize background quantities for the elliptic DG operator, possibly including the metric necessar...
Definition: Initialization.hpp:105
Tags::Normalized
Definition: Magnitude.hpp:137
Parallel::AlgorithmExecution::Continue
@ Continue
Leave the algorithm termination flag in its current state.
DataBox.hpp
elliptic::dg::Tags::MortarData
Holds elliptic::dg::MortarData, i.e. boundary data on both sides of a mortar.
Definition: DgOperator.hpp:223
elliptic::dg::Actions::apply_operator
tmpl::list< detail::PrepareAndSendMortarData< System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag >, detail::ReceiveMortarDataAndApplyOperator< System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag > > apply_operator
Apply the DG operator to the PrimalFieldsTag and write the result to the OperatorAppliedToFieldsTag
Definition: ApplyOperator.hpp:488
dg::SimpleBoundaryData
Distinguishes between field data, which can be projected to a mortar, and extra data,...
Definition: SimpleBoundaryData.hpp:32
cstddef
Assert.hpp
domain::Tags::InverseJacobian
The inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:165
elliptic::get_sources_computer
tmpl::conditional_t< Linearized, typename detail::sources_computer_linearized< System >::type, typename System::sources_computer > get_sources_computer
The System::sources_computer or the System::sources_computer_linearized, depending on the Linearized ...
Definition: GetSourcesComputer.hpp:30
elliptic::dg::Tags::PenaltyParameter
The prefactor to the penalty term of the numerical flux.
Definition: Tags.hpp:54
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::map
elliptic::dg::Actions::initialize_operator
tmpl::list< detail::InitializeFacesMortarsAndBackground< System, BackgroundTag > > initialize_operator
Initialize geometric and background quantities for the elliptic DG operator.
Definition: ApplyOperator.hpp:450
Parallel::AlgorithmExecution
AlgorithmExecution
The possible options for altering the current execution of the algorithm, used in the return type of ...
Definition: AlgorithmMetafunctions.hpp:31
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Element.hpp
DirectionMap
Definition: DirectionMap.hpp:15
ApplyAt.hpp
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
Tags::Magnitude
Definition: Magnitude.hpp:98
Gsl.hpp
db::mutate_apply
constexpr decltype(auto) mutate_apply(F &&f, const gsl::not_null< DataBox< BoxTags > * > box, Args &&... args) noexcept
Apply the invokable f mutating items MutateTags and taking as additional arguments ArgumentTags and a...
Definition: DataBox.hpp:1257
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
Frame::Logical
Definition: IndexType.hpp:42
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
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
elliptic::apply_boundary_condition
void apply_boundary_condition(const elliptic::BoundaryConditions::BoundaryCondition< Dim, Registrars > &boundary_condition, const db::DataBox< DbTagsList > &box, const MapKeys &map_keys_to_direction, const gsl::not_null< FieldsAndFluxes * >... fields_and_fluxes) noexcept
Apply the boundary_condition to the fields_and_fluxes with arguments from interface tags in the DataB...
Definition: ApplyBoundaryCondition.hpp:41
type_traits
TMPL.hpp
Mesh.hpp
Tags::Mortars
Definition: Tags.hpp:37
Parallel::InboxInserters::Map
Inserter for inserting data that is received as the value_type (with non-const key_type) of a map dat...
Definition: InboxInserters.hpp:36
string