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"
17 #include "Domain/InterfaceComputeTags.hpp"
18 #include "Domain/InterfaceHelpers.hpp"
19 #include "Domain/Structure/CreateInitialMesh.hpp"
22 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
23 #include "Domain/Tags.hpp"
24 #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
25 #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
26 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
27 #include "Elliptic/Systems/GetSourcesComputer.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
30 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
32 #include "Parallel/AlgorithmMetafunctions.hpp"
33 #include "Parallel/GlobalCache.hpp"
34 #include "Parallel/InboxInserters.hpp"
35 #include "Parallel/Invoke.hpp"
36 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
38 #include "Utilities/GetOutput.hpp"
39 #include "Utilities/Gsl.hpp"
40 #include "Utilities/Literals.hpp"
41 #include "Utilities/TMPL.hpp"
42 #include "Utilities/TaggedTuple.hpp"
43 
44 /// Actions related to elliptic discontinuous Galerkin schemes
46 // The individual actions in this namespace are not exposed publicly because
47 // they don't work on their own. Instead, the public interface (defined below)
48 // exposes them in action lists.
49 namespace detail {
50 
51 // This tag is used to communicate mortar data across neighbors
52 template <size_t Dim, typename TemporalIdTag, typename PrimalFields,
53  typename PrimalFluxes>
54 struct MortarDataInboxTag
56  MortarDataInboxTag<Dim, TemporalIdTag, PrimalFields, PrimalFluxes>> {
57  using temporal_id = typename TemporalIdTag::type;
58  using type = std::map<
59  temporal_id,
62  boost::hash<::dg::MortarId<Dim>>>>;
63 };
64 
65 // Initializes all quantities the DG operator needs on internal and external
66 // faces, as well as the mortars between neighboring elements.
67 template <typename System>
68 struct InitializeFacesAndMortars {
69  private:
70  static constexpr size_t Dim = System::volume_dim;
71 
72  // The DG operator needs these background fields on faces:
73  // - The `inv_metric_tag` on internal and external faces to normalize face
74  // normals (should be included in background fields)
75  // - The background fields in the `System::fluxes_computer::argument_tags` on
76  // internal and external faces to compute fluxes
77  // - All background fields on external faces to impose boundary conditions
78  using slice_tags_to_faces = tmpl::conditional_t<
79  std::is_same_v<typename System::background_fields, tmpl::list<>>,
80  tmpl::list<>,
81  tmpl::list<
82  // Possible optimization: Only the background fields in the
83  // System::fluxes_computer::argument_tags are needed on internal faces
85  using face_compute_tags = tmpl::list<
86  // For boundary conditions. Possible optimization: Not all systems need
87  // the coordinates on internal faces. Adding the compute item shouldn't
88  // incur a runtime overhead though since it's lazily evaluated.
90 
91  template <typename TagToSlice, typename DirectionsTag>
92  struct make_slice_tag {
94  };
95 
96  template <typename ComputeTag, typename DirectionsTag>
97  struct make_face_compute_tag {
99  };
100 
101  template <typename DirectionsTag>
102  using face_tags = tmpl::flatten<tmpl::list<
103  domain::Tags::InterfaceCompute<DirectionsTag,
105  domain::Tags::InterfaceCompute<DirectionsTag,
107  tmpl::transform<slice_tags_to_faces,
108  make_slice_tag<tmpl::_1, tmpl::pin<DirectionsTag>>>,
112  DirectionsTag,
113  tmpl::conditional_t<
114  std::is_same_v<typename System::inv_metric_tag, void>,
119  typename System::inv_metric_tag>>>,
121  DirectionsTag,
123  tmpl::transform<
125  make_face_compute_tag<tmpl::_1, tmpl::pin<DirectionsTag>>>>>;
126 
127  public:
128  using simple_tags =
129  tmpl::list<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
130  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>;
131  using compute_tags = tmpl::append<
132  tmpl::list<domain::Tags::InternalDirectionsCompute<Dim>,
134  face_tags<domain::Tags::InternalDirections<Dim>>,
135  face_tags<domain::Tags::BoundaryDirectionsInterior<Dim>>>;
136 
137  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
138  typename ArrayIndex, typename ActionList,
139  typename ParallelComponent>
140  static std::tuple<db::DataBox<DbTagsList>&&> apply(
141  db::DataBox<DbTagsList>& box,
142  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
143  const Parallel::GlobalCache<Metavariables>& /*cache*/,
144  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
145  const ParallelComponent* const /*meta*/) noexcept {
146  const auto& element = db::get<domain::Tags::Element<Dim>>(box);
147  const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
148  const auto& initial_extents =
149  db::get<domain::Tags::InitialExtents<Dim>>(box);
150 
151  // Setup initial mortars at interfaces to neighbors
152  ::dg::MortarMap<Dim, Mesh<Dim - 1>> all_mortar_meshes{};
153  ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>> all_mortar_sizes{};
154  for (const auto& [direction, neighbors] : element.neighbors()) {
155  const auto face_mesh = mesh.slice_away(direction.dimension());
156  const auto& orientation = neighbors.orientation();
157  for (const auto& neighbor : neighbors) {
158  const ::dg::MortarId<Dim> mortar_id{direction, neighbor};
159  all_mortar_meshes.emplace(
160  mortar_id,
161  ::dg::mortar_mesh(
162  face_mesh,
164  initial_extents, neighbor, mesh.quadrature(0), orientation)
165  .slice_away(direction.dimension())));
166  all_mortar_sizes.emplace(
167  mortar_id, ::dg::mortar_size(element.id(), neighbor,
168  direction.dimension(), orientation));
169  }
170  }
171  ::Initialization::mutate_assign<simple_tags>(make_not_null(&box),
172  std::move(all_mortar_meshes),
173  std::move(all_mortar_sizes));
174  return {std::move(box)};
175  }
176 };
177 
178 // Compute auxiliary variables and fluxes from the primal variables, prepare the
179 // local side of all mortars and send the local mortar data to neighbors. Also
180 // handle boundary conditions by preparing the exterior ("ghost") side of
181 // external mortars.
182 template <typename System, bool Linearized, typename TemporalIdTag,
183  typename PrimalFieldsTag, typename PrimalFluxesTag,
184  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
185  typename PrimalMortarFluxesTag,
186  typename FluxesArgsTags =
187  typename System::fluxes_computer::argument_tags,
188  typename SourcesArgsTags = typename elliptic::get_sources_computer<
189  System, Linearized>::argument_tags>
190 struct PrepareAndSendMortarData;
191 
192 template <typename System, bool Linearized, typename TemporalIdTag,
193  typename PrimalFieldsTag, typename PrimalFluxesTag,
194  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
195  typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
196  typename... SourcesArgsTags>
197 struct PrepareAndSendMortarData<
198  System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
199  OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
200  tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
201  private:
202  static constexpr size_t Dim = System::volume_dim;
203  using all_mortar_data_tag = ::Tags::Mortars<
204  elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
205  typename PrimalMortarFieldsTag::tags_list,
206  typename PrimalMortarFluxesTag::tags_list>,
207  Dim>;
208  using mortar_data_inbox_tag =
209  MortarDataInboxTag<Dim, TemporalIdTag,
210  typename PrimalMortarFieldsTag::tags_list,
211  typename PrimalMortarFluxesTag::tags_list>;
212  using BoundaryConditionsBase = typename System::boundary_conditions_base;
213 
214  public:
215  // Request these tags be added to the DataBox by the `SetupDataBox` action. We
216  // don't actually need to initialize them, because the `TemporalIdTag` and the
217  // `PrimalFieldsTag` will be set by other actions before applying the operator
218  // and the remaining tags hold output of the operator.
219  using simple_tags =
220  tmpl::list<TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
221  OperatorAppliedToFieldsTag, all_mortar_data_tag>;
222  using compute_tags = tmpl::list<>;
223 
224  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
225  typename ActionList, typename ParallelComponent>
226  static std::tuple<db::DataBox<DbTagsList>&&> apply(
227  db::DataBox<DbTagsList>& box,
228  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
230  const ElementId<Dim>& element_id, const ActionList /*meta*/,
231  const ParallelComponent* const /*meta*/) noexcept {
232  const auto& temporal_id = db::get<TemporalIdTag>(box);
233  const auto& element = db::get<domain::Tags::Element<Dim>>(box);
234  const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
235  const size_t num_points = mesh.number_of_grid_points();
236  const auto& mortar_meshes =
237  db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box);
238  const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
239  const auto& boundary_conditions = domain.blocks()
240  .at(element_id.block_id())
241  .external_boundary_conditions();
242  const auto apply_boundary_condition =
243  [&box, &boundary_conditions, &element_id](
244  const Direction<Dim>& direction,
245  const auto... fields_and_fluxes) noexcept {
246  ASSERT(
247  boundary_conditions.contains(direction),
248  "No boundary condition is available in block "
249  << element_id.block_id() << " in direction " << direction
250  << ". Make sure you are setting up boundary conditions when "
251  "creating the domain.");
252  ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
253  boundary_conditions.at(direction).get()) != nullptr,
254  "The boundary condition in block "
255  << element_id.block_id() << " in direction " << direction
256  << " has an unexpected type. Make sure it derives off the "
257  "'boundary_conditions_base' class set in the system.");
258  const auto& boundary_condition =
259  dynamic_cast<const BoundaryConditionsBase&>(
260  *boundary_conditions.at(direction));
261  elliptic::apply_boundary_condition<Linearized>(
262  boundary_condition, box, direction, fields_and_fluxes...);
263  };
264 
265  // Can't `db::get` the arguments for the boundary conditions within
266  // `db::mutate`, so we extract the data to mutate and move it back in when
267  // we're done.
268  typename PrimalFluxesTag::type primal_fluxes;
269  typename all_mortar_data_tag::type all_mortar_data;
270  db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
271  make_not_null(&box),
272  [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
273  const auto local_all_mortar_data) {
274  primal_fluxes = std::move(*local_primal_fluxes);
275  all_mortar_data = std::move(*local_all_mortar_data);
276  });
277 
278  // Prepare mortar data
279  //
280  // These memory buffers will be discarded when the action returns so we
281  // don't inflate the memory usage of the simulation when the element is
282  // inactive.
283  Variables<typename System::auxiliary_fields> auxiliary_fields_buffer{
284  num_points};
285  Variables<typename System::auxiliary_fluxes> auxiliary_fluxes_buffer{
286  num_points};
287  elliptic::dg::prepare_mortar_data<System, Linearized>(
288  make_not_null(&auxiliary_fields_buffer),
289  make_not_null(&auxiliary_fluxes_buffer), make_not_null(&primal_fluxes),
290  make_not_null(&all_mortar_data), db::get<PrimalFieldsTag>(box), element,
293  Frame::Inertial>>(box),
297  box),
301  box),
308  mortar_meshes,
310  temporal_id, apply_boundary_condition,
311  std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
312  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
314  typename System::fluxes_computer::argument_tags,
315  get_volume_tags<typename System::fluxes_computer>>(
316  [](const auto&... fluxes_args_on_face) noexcept {
317  return std::forward_as_tuple(fluxes_args_on_face...);
318  },
319  box),
321  typename System::fluxes_computer::argument_tags,
322  get_volume_tags<typename System::fluxes_computer>>(
323  [](const auto&... fluxes_args_on_face) noexcept {
324  return std::forward_as_tuple(fluxes_args_on_face...);
325  },
326  box));
327 
328  // Move the mutated data back into the DataBox
329  db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
330  make_not_null(&box),
331  [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
332  const auto local_all_mortar_data) {
333  *local_primal_fluxes = std::move(primal_fluxes);
334  *local_all_mortar_data = std::move(all_mortar_data);
335  });
336 
337  // Send mortar data to neighbors
338  auto& receiver_proxy =
339  Parallel::get_parallel_component<ParallelComponent>(cache);
340  for (const auto& [direction, neighbors] : element.neighbors()) {
341  const size_t dimension = direction.dimension();
342  const auto& orientation = neighbors.orientation();
343  const auto direction_from_neighbor = orientation(direction.opposite());
344  for (const auto& neighbor_id : neighbors) {
345  const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
346  // Make a copy of the local boundary data on the mortar to send to the
347  // neighbor
348  auto remote_boundary_data_on_mortar =
349  get<all_mortar_data_tag>(box).at(mortar_id).local_data(
350  temporal_id);
351  // Reorient the data to the neighbor orientation if necessary
352  if (not orientation.is_aligned()) {
353  remote_boundary_data_on_mortar.orient_on_slice(
354  mortar_meshes.at(mortar_id).extents(), dimension, orientation);
355  }
356  // Send remote data to neighbor
357  Parallel::receive_data<mortar_data_inbox_tag>(
358  receiver_proxy[neighbor_id], temporal_id,
359  std::make_pair(
360  ::dg::MortarId<Dim>{direction_from_neighbor, element.id()},
361  std::move(remote_boundary_data_on_mortar)));
362  } // loop over neighbors in direction
363  } // loop over directions
364 
365  return {std::move(box)};
366  }
367 };
368 
369 // Wait until all mortar data from neighbors is available. Then add boundary
370 // corrections to the primal fluxes, compute their derivatives (i.e. the second
371 // derivatives of the primal variables) and add boundary corrections to the
372 // result.
373 template <typename System, bool Linearized, typename TemporalIdTag,
374  typename PrimalFieldsTag, typename PrimalFluxesTag,
375  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
376  typename PrimalMortarFluxesTag,
377  typename FluxesArgsTags =
378  typename System::fluxes_computer::argument_tags,
379  typename SourcesArgsTags = typename elliptic::get_sources_computer<
380  System, Linearized>::argument_tags>
381 struct ReceiveMortarDataAndApplyOperator;
382 
383 template <typename System, bool Linearized, typename TemporalIdTag,
384  typename PrimalFieldsTag, typename PrimalFluxesTag,
385  typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
386  typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
387  typename... SourcesArgsTags>
388 struct ReceiveMortarDataAndApplyOperator<
389  System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
390  OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
391  tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
392  private:
393  static constexpr size_t Dim = System::volume_dim;
394  using all_mortar_data_tag = ::Tags::Mortars<
395  elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
396  typename PrimalMortarFieldsTag::tags_list,
397  typename PrimalMortarFluxesTag::tags_list>,
398  Dim>;
399  using mortar_data_inbox_tag =
400  MortarDataInboxTag<Dim, TemporalIdTag,
401  typename PrimalMortarFieldsTag::tags_list,
402  typename PrimalMortarFluxesTag::tags_list>;
403 
404  public:
405  using const_global_cache_tags =
406  tmpl::list<elliptic::dg::Tags::PenaltyParameter>;
407  using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
408 
409  template <typename DbTags, typename... InboxTags, typename Metavariables,
410  typename ArrayIndex, typename ActionList,
411  typename ParallelComponent>
413  db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
414  const Parallel::GlobalCache<Metavariables>& /*cache*/,
415  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
416  const ParallelComponent* const /*meta*/) noexcept {
417  const auto& temporal_id = get<TemporalIdTag>(box);
418 
419  if (not ::dg::has_received_from_all_mortars<mortar_data_inbox_tag>(
420  temporal_id, get<domain::Tags::Element<Dim>>(box), inboxes)) {
421  return {std::move(box), Parallel::AlgorithmExecution::Retry};
422  }
423 
424  // Move received "remote" mortar data into the DataBox
425  if (LIKELY(db::get<domain::Tags::Element<Dim>>(box).number_of_neighbors() >
426  0)) {
427  auto received_mortar_data =
428  std::move(tuples::get<mortar_data_inbox_tag>(inboxes)
429  .extract(temporal_id)
430  .mapped());
431  db::mutate<all_mortar_data_tag>(
432  make_not_null(&box), [&received_mortar_data, &temporal_id](
433  const auto all_mortar_data) noexcept {
434  for (auto& [mortar_id, mortar_data] : received_mortar_data) {
435  all_mortar_data->at(mortar_id).remote_insert(
436  temporal_id, std::move(mortar_data));
437  }
438  });
439  }
440 
441  // Apply DG operator
442  db::mutate<OperatorAppliedToFieldsTag, all_mortar_data_tag>(
443  make_not_null(&box),
444  [](const auto&... args) noexcept {
445  elliptic::dg::apply_operator<System, Linearized>(args...);
446  },
447  db::get<PrimalFieldsTag>(box), db::get<PrimalFluxesTag>(box),
448  db::get<domain::Tags::Mesh<Dim>>(box),
450  Frame::Inertial>>(box),
457  db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
459  db::get<elliptic::dg::Tags::PenaltyParameter>(box), temporal_id,
460  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
461 
462  return {std::move(box), Parallel::AlgorithmExecution::Continue};
463  }
464 };
465 
466 } // namespace detail
467 
468 /*!
469  * \brief Initialize DataBox tags for applying the DG operator
470  *
471  * \see elliptic::dg::Actions::apply_operator
472  */
473 template <typename System>
474 using initialize_operator =
475  tmpl::list<detail::InitializeFacesAndMortars<System>>;
476 
477 /*!
478  * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
479  * the `OperatorAppliedToFieldsTag`
480  *
481  * Add this list to the action list of a parallel component to compute the
482  * elliptic DG operator or its linearization. The operator involves a
483  * communication between nearest-neighbor elements. See `elliptic::dg` for
484  * details on the elliptic DG operator. Make sure to add
485  * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
486  * your parallel component so the required DataBox tags are set up before
487  * applying the operator.
488  *
489  * The result of the computation is written to the `OperatorAppliedToFieldsTag`.
490  * Additionally, the primal fluxes are written to the `PrimalFluxesTag` as an
491  * intermediate result. The auxiliary fields and fluxes are discarded to avoid
492  * inflating the memory usage.
493  *
494  * You can specify the `PrimalMortarFieldsTag` and the `PrimalMortarFluxesTag`
495  * to re-use mortar-data memory buffers from other operator applications, for
496  * example when applying the nonlinear and linearized operator. They default to
497  * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
498  * corresponding to these tags are set up in the DataBox.
499  */
500 template <typename System, bool Linearized, typename TemporalIdTag,
501  typename PrimalFieldsTag, typename PrimalFluxesTag,
502  typename OperatorAppliedToFieldsTag,
503  typename PrimalMortarFieldsTag = PrimalFieldsTag,
504  typename PrimalMortarFluxesTag = PrimalFluxesTag>
505 using apply_operator =
506  tmpl::list<detail::PrepareAndSendMortarData<
507  System, Linearized, TemporalIdTag, PrimalFieldsTag,
508  PrimalFluxesTag, OperatorAppliedToFieldsTag,
509  PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
510  detail::ReceiveMortarDataAndApplyOperator<
511  System, Linearized, TemporalIdTag, PrimalFieldsTag,
512  PrimalFluxesTag, OperatorAppliedToFieldsTag,
513  PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
514 
515 /*!
516  * \brief For linear systems, impose inhomogeneous boundary conditions as
517  * contributions to the fixed sources (i.e. the RHS of the equations).
518  *
519  * \see elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
520  */
521 template <
522  typename System, typename FixedSourcesTag,
523  typename FluxesArgsTags = typename System::fluxes_computer::argument_tags,
524  typename SourcesArgsTags = typename System::sources_computer::argument_tags>
526 
527 /// \cond
528 template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
529  typename... SourcesArgsTags>
531  System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
532  tmpl::list<SourcesArgsTags...>> {
533  private:
534  static constexpr size_t Dim = System::volume_dim;
535  using BoundaryConditionsBase = typename System::boundary_conditions_base;
536 
537  public:
538  using const_global_cache_tags =
539  tmpl::list<elliptic::dg::Tags::PenaltyParameter>;
540 
541  template <typename DbTags, typename... InboxTags, typename Metavariables,
542  typename ActionList, typename ParallelComponent>
543  static std::tuple<db::DataBox<DbTags>&&> apply(
544  db::DataBox<DbTags>& box,
545  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
546  const Parallel::GlobalCache<Metavariables>& /*cache*/,
547  const ElementId<Dim>& element_id, const ActionList /*meta*/,
548  const ParallelComponent* const /*meta*/) noexcept {
549  const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
550  const auto& boundary_conditions = domain.blocks()
551  .at(element_id.block_id())
552  .external_boundary_conditions();
553  const auto apply_boundary_condition =
554  [&box, &boundary_conditions, &element_id](
555  const Direction<Dim>& direction,
556  const auto... fields_and_fluxes) noexcept {
557  ASSERT(
558  boundary_conditions.contains(direction),
559  "No boundary condition is available in block "
560  << element_id.block_id() << " in direction " << direction
561  << ". Make sure you are setting up boundary conditions when "
562  "creating the domain.");
563  ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
564  boundary_conditions.at(direction).get()) != nullptr,
565  "The boundary condition in block "
566  << element_id.block_id() << " in direction " << direction
567  << " has an unexpected type. Make sure it derives off the "
568  "'boundary_conditions_base' class set in the system.");
569  const auto& boundary_condition =
570  dynamic_cast<const BoundaryConditionsBase&>(
571  *boundary_conditions.at(direction));
572  elliptic::apply_boundary_condition<false>(
573  boundary_condition, box, direction, fields_and_fluxes...);
574  };
575 
576  // Can't `db::get` the arguments for the boundary conditions within
577  // `db::mutate`, so we extract the data to mutate and move it back in when
578  // we're done.
579  typename FixedSourcesTag::type fixed_sources;
580  db::mutate<FixedSourcesTag>(
581  make_not_null(&box), [&fixed_sources](const auto local_fixed_sources) {
582  fixed_sources = std::move(*local_fixed_sources);
583  });
584 
585  elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source<System>(
586  make_not_null(&fixed_sources), db::get<domain::Tags::Element<Dim>>(box),
589  Frame::Inertial>>(box),
593  box),
599  db::get<elliptic::dg::Tags::PenaltyParameter>(box),
601  std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
602  std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
604  typename System::fluxes_computer::argument_tags,
605  get_volume_tags<typename System::fluxes_computer>>(
606  [](const auto&... fluxes_args_on_face) noexcept {
607  return std::forward_as_tuple(fluxes_args_on_face...);
608  },
609  box),
611  typename System::fluxes_computer::argument_tags,
612  get_volume_tags<typename System::fluxes_computer>>(
613  [](const auto&... fluxes_args_on_face) noexcept {
614  return std::forward_as_tuple(fluxes_args_on_face...);
615  },
616  box));
617 
618  // Move the mutated data back into the DataBox
619  db::mutate<FixedSourcesTag>(
620  make_not_null(&box), [&fixed_sources](const auto local_fixed_sources) {
621  *local_fixed_sources = std::move(fixed_sources);
622  });
623  return {std::move(box)};
624  }
625 };
626 /// \endcond
627 
628 } // namespace elliptic::dg::Actions
domain::Tags::UnnormalizedFaceNormal
Definition: FaceNormal.hpp:112
elliptic::dg::Actions
Actions related to elliptic discontinuous Galerkin schemes.
Definition: ApplyOperator.hpp:45
maximum_number_of_neighbors
constexpr size_t maximum_number_of_neighbors(const size_t dim)
Definition: MaxNumberOfNeighbors.hpp:15
elliptic::apply_boundary_condition
void apply_boundary_condition(const elliptic::BoundaryConditions::BoundaryCondition< Dim, Registrars > &boundary_condition, const db::DataBox< DbTagsList > &box, const Direction< Dim > &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:36
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 ...
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
Literals.hpp
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: ApplyOperator.hpp:525
std::pair
GlobalCache.hpp
Tags::MortarSize
Definition: Tags.hpp:48
FaceNormal.hpp
Tags.hpp
domain::Tags::Element
Definition: Tags.hpp:97
domain::Tags::UnnormalizedFaceNormalCompute
Definition: FaceNormal.hpp:117
Tags::Variables
Definition: VariablesTag.hpp:21
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
dg::Initialization::face_compute_tags
tmpl::list< Tags... > face_compute_tags
Definition: InitializeInterfaces.hpp:41
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:297
Direction< Dim >
domain::Tags::InterfaceMesh
Definition: InterfaceComputeTags.hpp:167
ElementId< Dim >
domain::Tags::Slice
Compute tag for representing a compute item that slices data from the volume to a set of interfaces.
Definition: InterfaceComputeTags.hpp:115
Tags::Normalized
Definition: Magnitude.hpp:137
Tags::EuclideanMagnitude
Definition: Magnitude.hpp:109
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:513
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
std::array
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::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
Tags::NonEuclideanMagnitude
Definition: Magnitude.hpp:124
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
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
domain::Tags::Interface
Tag which is either a SimpleTag for quantities on an interface, base tag to a compute item which acts...
Definition: Tags.hpp:366
domain::Initialization::create_initial_mesh
Mesh< Dim > create_initial_mesh(const std::vector< std::array< size_t, Dim >> &initial_extents, const ElementId< Dim > &element_id, Spectral::Quadrature quadrature, const OrientationMap< Dim > &orientation={}) noexcept
Construct the initial Mesh of an Element.
domain::Tags::InterfaceCompute
Compute tag for representing items computed on a set of interfaces. Can be retrieved using Tags::Inte...
Definition: InterfaceComputeTags.hpp:67
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
domain::Tags::Direction
Definition: Tags.hpp:382
domain::Tags::InternalDirections
Definition: Tags.hpp:270
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
Frame::Logical
Definition: IndexType.hpp:42
Tags::NormalizedCompute
Definition: Magnitude.hpp:148
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::BoundaryDirectionsInteriorCompute
Definition: Tags.hpp:303
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
Mesh.hpp
domain::Tags::BoundaryCoordinates
Definition: InterfaceComputeTags.hpp:186
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
elliptic::dg::Actions::initialize_operator
tmpl::list< detail::InitializeFacesAndMortars< System > > initialize_operator
Initialize DataBox tags for applying the DG operator.
Definition: ApplyOperator.hpp:475
string