SubdomainOperator.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <optional>
8 #include <ostream>
9 #include <tuple>
10 #include <utility>
11 
13 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
15 #include "Domain/Domain.hpp"
16 #include "Domain/FaceNormal.hpp"
17 #include "Domain/InterfaceHelpers.hpp"
19 #include "Domain/Structure/DirectionMap.hpp"
21 #include "Domain/Tags.hpp"
22 #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
23 #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
24 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
25 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
26 #include "Elliptic/Systems/GetSourcesComputer.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
31 #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
32 #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
33 #include "ParallelAlgorithms/LinearSolver/Schwarz/SubdomainOperator.hpp"
34 #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
36 #include "Utilities/Gsl.hpp"
37 #include "Utilities/TMPL.hpp"
38 
39 /// Items related to the restriction of the DG operator to an element-centered
40 /// subdomain
42 
43 namespace detail {
44 // Wrap the `Tag` in `LinearSolver::Schwarz::Tags::Overlaps`, except if it is
45 // included in `TakeFromCenterTags`.
46 template <typename Tag, typename Dim, typename OptionsGroup,
47  typename TakeFromCenterTags>
48 struct make_overlap_tag_impl {
49  using type = tmpl::conditional_t<
50  tmpl::list_contains_v<TakeFromCenterTags, Tag>, Tag,
52 };
53 
54 // Wrap the `Tag` in `Tags::NeighborMortars`
55 template <typename Tag, typename Dim>
56 struct make_neighbor_mortars_tag_impl {
57  using type = Tags::NeighborMortars<Tag, Dim::value>;
58 };
59 } // namespace detail
60 
61 /*!
62  * \brief The elliptic DG operator on an element-centered subdomain
63  *
64  * This operator is a restriction of the full (linearized) DG-operator to an
65  * element-centered subdomain with a few points overlap into neighboring
66  * elements. It is a `LinearSolver::Schwarz::SubdomainOperator` to be used with
67  * the Schwarz linear solver when it solves the elliptic DG operator.
68  *
69  * This operator requires the following tags are available on overlap regions
70  * with neighboring elements:
71  *
72  * - Geometric quantities provided by
73  * `elliptic::dg::subdomain_operator::InitializeSubdomain`.
74  * - All `System::fluxes_computer::argument_tags` and
75  * `System::sources_computer::argument_tags` (or
76  * `System::sources_computer_linearized::argument_tags` for nonlinear
77  * systems), except those listed in `ArgsTagsFromCenter`. The latter will be
78  * taken from the central element's DataBox, so they don't need to be made
79  * available on overlaps.
80  * - The `System::fluxes_computer::argument_tags` on internal and external
81  * interfaces, except those listed in `System::fluxes_computer::volume_tags`.
82  *
83  * Some of these tags may require communication between elements. For example,
84  * nonlinear system fields are constant background fields for the linearized DG
85  * operator, but are updated in every nonlinear solver iteration. Therefore, the
86  * updated nonlinear fields must be communicated across overlaps between
87  * nonlinear solver iterations. To perform the communication you can use
88  * `LinearSolver::Schwarz::Actions::SendOverlapFields` and
89  * `LinearSolver::Schwarz::Actions::ReceiveOverlapFields`, setting
90  * `RestrictToOverlap` to `false`. See
91  * `LinearSolver::Schwarz::SubdomainOperator` for details.
92  *
93  * \warning The subdomain operator hasn't been tested with periodic boundary
94  * conditions so far.
95  */
96 template <typename System, typename OptionsGroup,
97  typename ArgsTagsFromCenter = tmpl::list<>>
99  : LinearSolver::Schwarz::SubdomainOperator<System::volume_dim> {
100  private:
101  static constexpr size_t Dim = System::volume_dim;
102 
103  // Operator applications happen sequentially so we don't have to keep track of
104  // the temporal id
105  static constexpr size_t temporal_id = 0;
106 
107  // The subdomain operator always applies the linearized DG operator
108  static constexpr bool linearized = true;
109 
110  using BoundaryConditionsBase = typename System::boundary_conditions_base;
111 
112  // These are the arguments that we need to retrieve from the DataBox and pass
113  // to the functions in `elliptic::dg`, both on the central element and on
114  // neighbors
115  using prepare_args_tags = tmpl::list<
130  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
131  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>;
132  using apply_args_tags = tmpl::list<
142  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
143  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
145  using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
146  using sources_args_tags =
147  typename elliptic::get_sources_computer<System,
148  linearized>::argument_tags;
149 
150  // We need the fluxes args also on interfaces (internal and external). The
151  // volume tags are the subset that don't have to be taken from interfaces.
152  using fluxes_args_volume_tags =
153  get_volume_tags<typename System::fluxes_computer>;
154 
155  // These tags can be taken directly from the central element's DataBox, even
156  // when evaluating neighbors
157  using args_tags_from_center = tmpl::remove_duplicates<
160 
161  // Data on neighbors is stored in the central element's DataBox in
162  // `LinearSolver::Schwarz::Tags::Overlaps` maps, so we wrap the argument tags
163  // with this prefix
164  using make_overlap_tag =
165  detail::make_overlap_tag_impl<tmpl::_1, tmpl::pin<tmpl::size_t<Dim>>,
166  tmpl::pin<OptionsGroup>,
167  tmpl::pin<args_tags_from_center>>;
168  using prepare_args_tags_overlap =
169  tmpl::transform<prepare_args_tags, make_overlap_tag>;
170  using apply_args_tags_overlap =
171  tmpl::transform<apply_args_tags, make_overlap_tag>;
172  using fluxes_args_tags_overlap =
173  tmpl::transform<fluxes_args_tags, make_overlap_tag>;
174  using sources_args_tags_overlap =
175  tmpl::transform<sources_args_tags, make_overlap_tag>;
176  template <typename Directions>
177  using fluxes_args_tags_overlap_interface = tmpl::transform<
178  tmpl::transform<fluxes_args_tags,
180  tmpl::pin<fluxes_args_volume_tags>>>,
181  make_overlap_tag>;
182 
183  // We also need some data on the remote side of all neighbors' mortars. Such
184  // data is stored in the central element's DataBox in `Tags::NeighborMortars`
185  // maps
186  using make_neighbor_mortars_tag =
187  detail::make_neighbor_mortars_tag_impl<tmpl::_1,
188  tmpl::pin<tmpl::size_t<Dim>>>;
189 
190  public:
191  template <typename ResultTags, typename OperandTags, typename DbTagsList>
192  void operator()(
193  const gsl::not_null<
195  result,
197  Dim, OperandTags>& operand,
198  db::DataBox<DbTagsList>& box) noexcept {
199  // Used to retrieve items out of the DataBox to forward to functions. This
200  // replaces a long series of db::get calls.
201  const auto get_items = [](const auto&... args) noexcept {
202  return std::forward_as_tuple(args...);
203  };
204 
205  // Retrieve data out of the DataBox
206  using tags_to_retrieve = tmpl::flatten<tmpl::list<
208  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
209  // Data on overlaps with neighbors
210  tmpl::transform<
211  tmpl::flatten<tmpl::list<
218  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
219  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
220  // Data on the remote side of the neighbor's mortars
221  tmpl::transform<
222  tmpl::list<domain::Tags::Mesh<Dim>,
225  domain::Tags::Mesh<Dim - 1>,
226  ::Tags::MortarSize<Dim - 1>>,
227  make_neighbor_mortars_tag>>>,
228  make_overlap_tag>>>;
229  const auto& [domain, central_element, central_mortar_meshes,
230  all_overlap_extents, all_neighbor_elements,
231  all_neighbor_meshes,
232  all_neighbor_face_normal_magnitudes_internal,
233  all_neighbor_mortar_meshes, all_neighbor_mortar_sizes,
234  all_neighbors_neighbor_meshes,
235  all_neighbors_neighbor_face_normal_magnitudes,
236  all_neighbors_neighbor_mortar_meshes,
237  all_neighbors_neighbor_mortar_sizes] =
238  db::apply<tags_to_retrieve>(get_items, box);
239  const auto fluxes_args = db::apply<fluxes_args_tags>(get_items, box);
240  const auto sources_args = db::apply<sources_args_tags>(get_items, box);
241  const auto fluxes_args_on_internal_faces =
242  interface_apply<domain::Tags::InternalDirections<Dim>, fluxes_args_tags,
243  fluxes_args_volume_tags>(get_items, box);
244  const auto fluxes_args_on_external_faces =
245  interface_apply<domain::Tags::BoundaryDirectionsInterior<Dim>,
246  fluxes_args_tags, fluxes_args_volume_tags>(get_items,
247  box);
248 
249  // Setup boundary conditions
250  const auto apply_boundary_condition =
251  [&box, &local_domain = domain](
252  const ElementId<Dim>& local_element_id,
253  const Direction<Dim>& local_direction, auto is_overlap,
254  const auto& map_keys, const auto... fields_and_fluxes) noexcept {
255  constexpr bool is_overlap_v =
256  std::decay_t<decltype(is_overlap)>::value;
257  const auto& boundary_conditions = local_domain.blocks()
258  .at(local_element_id.block_id())
259  .external_boundary_conditions();
260  ASSERT(boundary_conditions.contains(local_direction),
261  "No boundary condition is available in block "
262  << local_element_id.block_id() << " in direction "
263  << local_direction
264  << ". Make sure you are setting up boundary conditions "
265  "when creating the domain.");
266  ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
267  boundary_conditions.at(local_direction).get()) != nullptr,
268  "The boundary condition in block "
269  << local_element_id.block_id() << " in direction "
270  << local_direction
271  << " has an unexpected type. Make sure it derives off the "
272  "'boundary_conditions_base' class set in the system.");
273  const auto& boundary_condition =
274  dynamic_cast<const BoundaryConditionsBase&>(
275  *boundary_conditions.at(local_direction));
277  linearized,
278  tmpl::conditional_t<is_overlap_v, make_overlap_tag, void>>(
279  boundary_condition, box, map_keys, fields_and_fluxes...);
280  };
281 
282  // The subdomain operator essentially does two sweeps over all elements in
283  // the subdomain: In the first sweep it prepares the mortar data and stores
284  // them on both sides of all mortars, and in the second sweep it consumes
285  // the mortar data to apply the operator. This implementation is relatively
286  // simple because it can re-use the implementation for the parallel DG
287  // operator. However, it is also possible to apply the subdomain operator in
288  // a single sweep over all elements, incrementally building up the mortar
289  // data and applying boundary corrections immediately to both adjacent
290  // elements once the data is available. That approach is possibly a
291  // performance optimization but requires re-implementing a lot of logic for
292  // the DG operator here. It should be considered once the subdomain operator
293  // has been identified as the performance bottleneck. An alternative to
294  // optimizing the subdomain operator performance is to precondition the
295  // subdomain solve with a _much_ simpler subdomain operator, such as a
296  // finite-difference Laplacian, so fewer applications of the more expensive
297  // DG subdomain operator are necessary.
298 
299  // 1. Prepare mortar data on all elements in the subdomain and store them on
300  // mortars, reorienting if needed
301  //
302  // Prepare central element
303  const auto apply_boundary_condition_center =
304  [&apply_boundary_condition, &local_central_element = central_element](
305  const Direction<Dim>& local_direction,
306  const auto... fields_and_fluxes) noexcept {
307  apply_boundary_condition(local_central_element.id(), local_direction,
308  std::false_type{}, local_direction,
309  fields_and_fluxes...);
310  };
311  db::apply<prepare_args_tags>(
312  [this, &operand](const auto&... args) noexcept {
313  elliptic::dg::prepare_mortar_data<System, linearized>(
314  make_not_null(&central_auxiliary_vars_),
315  make_not_null(&central_auxiliary_fluxes_),
316  make_not_null(&central_primal_fluxes_),
317  make_not_null(&central_mortar_data_), operand.element_data,
318  args...);
319  },
320  box, temporal_id, apply_boundary_condition_center, fluxes_args,
321  sources_args, fluxes_args_on_internal_faces,
322  fluxes_args_on_external_faces);
323  // Prepare neighbors
324  for (const auto& [direction, neighbors] : central_element.neighbors()) {
325  const auto& orientation = neighbors.orientation();
326  const auto direction_from_neighbor = orientation(direction.opposite());
327  for (const auto& neighbor_id : neighbors) {
328  const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
329  neighbor_id};
330  const auto& overlap_extent = all_overlap_extents.at(overlap_id);
331  const auto& neighbor = all_neighbor_elements.at(overlap_id);
332  const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
333  const auto& mortar_id = overlap_id;
334  const auto& mortar_mesh = central_mortar_meshes.at(mortar_id);
335  const ::dg::MortarId<Dim> mortar_id_from_neighbor{
336  direction_from_neighbor, central_element.id()};
337 
338  // Intercept empty overlaps. In the unlikely case that overlaps have
339  // zero extent, meaning no point of the neighbor is part of the
340  // subdomain (which is fairly useless, except for testing), the
341  // subdomain is identical to the central element and no communication
342  // with neighbors is necessary. We can just handle the mortar between
343  // central element and neighbor and continue.
344  if (UNLIKELY(overlap_extent == 0)) {
345  const auto& mortar_mesh_from_neighbor =
346  all_neighbor_mortar_meshes.at(overlap_id)
347  .at(mortar_id_from_neighbor);
348  const auto& mortar_size_from_neighbor =
349  all_neighbor_mortar_sizes.at(overlap_id)
350  .at(mortar_id_from_neighbor);
351  auto remote_boundary_data =
353  typename System::primal_fields,
354  typename System::primal_fluxes>(
355  direction_from_neighbor, neighbor_mesh,
356  all_neighbor_face_normal_magnitudes_internal.at(overlap_id)
357  .at(direction_from_neighbor),
358  mortar_mesh_from_neighbor, mortar_size_from_neighbor);
359  if (not orientation.is_aligned()) {
360  remote_boundary_data.orient_on_slice(
361  mortar_mesh_from_neighbor.extents(),
362  direction_from_neighbor.dimension(), orientation.inverse_map());
363  }
364  central_mortar_data_.at(mortar_id).remote_insert(
365  temporal_id, std::move(remote_boundary_data));
366  continue;
367  }
368 
369  // Copy the central element's mortar data to the neighbor
370  auto oriented_mortar_data =
371  central_mortar_data_.at(mortar_id).local_data(temporal_id);
372  if (not orientation.is_aligned()) {
373  oriented_mortar_data.orient_on_slice(
374  mortar_mesh.extents(), direction.dimension(), orientation);
375  }
376  neighbors_mortar_data_[overlap_id][::dg::MortarId<Dim>{
377  direction_from_neighbor,
378  central_element.id()}]
379  .remote_insert(temporal_id, std::move(oriented_mortar_data));
380 
381  // Now we switch perspective to the neighbor. First, we extend the
382  // overlap data to the full neighbor mesh by padding it with zeros. This
383  // is necessary because spectral operators such as derivatives require
384  // data on the full mesh.
386  make_not_null(&extended_operand_vars_[overlap_id]),
387  operand.overlap_data.at(overlap_id), neighbor_mesh.extents(),
388  overlap_extent, direction_from_neighbor);
389 
390  const auto apply_boundary_condition_neighbor =
391  [&apply_boundary_condition, &local_neighbor_id = neighbor_id,
392  &overlap_id](const Direction<Dim>& local_direction,
393  const auto... fields_and_fluxes) noexcept {
395  local_neighbor_id, local_direction, std::true_type{},
396  std::forward_as_tuple(overlap_id, local_direction),
397  fields_and_fluxes...);
398  };
399 
400  using FluxesArgs = std::decay_t<decltype(fluxes_args)>;
401  const auto fluxes_args_on_overlap =
402  elliptic::util::apply_at<fluxes_args_tags_overlap,
403  args_tags_from_center>(get_items, box,
404  overlap_id);
405  const auto sources_args_on_overlap =
406  elliptic::util::apply_at<sources_args_tags_overlap,
407  args_tags_from_center>(get_items, box,
408  overlap_id);
409  DirectionMap<Dim, FluxesArgs> fluxes_args_on_overlap_faces_internal{};
410  DirectionMap<Dim, FluxesArgs> fluxes_args_on_overlap_faces_external{};
411  for (const auto& neighbor_direction : neighbor.internal_boundaries()) {
412  fluxes_args_on_overlap_faces_internal.emplace(
413  neighbor_direction,
414  elliptic::util::apply_at<
415  fluxes_args_tags_overlap_interface<
417  args_tags_from_center>(
418  get_items, box,
419  std::forward_as_tuple(overlap_id, neighbor_direction)));
420  }
421  for (const auto& neighbor_direction : neighbor.external_boundaries()) {
422  fluxes_args_on_overlap_faces_external.emplace(
423  neighbor_direction,
424  elliptic::util::apply_at<
425  fluxes_args_tags_overlap_interface<
427  args_tags_from_center>(
428  get_items, box,
429  std::forward_as_tuple(overlap_id, neighbor_direction)));
430  }
431 
432  elliptic::util::apply_at<prepare_args_tags_overlap,
433  args_tags_from_center>(
434  [this, &overlap_id](const auto&... args) noexcept {
435  elliptic::dg::prepare_mortar_data<System, linearized>(
436  make_not_null(&neighbors_auxiliary_vars_[overlap_id]),
437  make_not_null(&neighbors_auxiliary_fluxes_[overlap_id]),
438  make_not_null(&neighbors_primal_fluxes_[overlap_id]),
439  make_not_null(&neighbors_mortar_data_[overlap_id]),
440  extended_operand_vars_[overlap_id], args...);
441  },
442  box, overlap_id, temporal_id, apply_boundary_condition_neighbor,
443  fluxes_args_on_overlap, sources_args_on_overlap,
444  fluxes_args_on_overlap_faces_internal,
445  fluxes_args_on_overlap_faces_external);
446 
447  // Copy this neighbor's mortar data to the other side of the mortars. On
448  // the other side we either have the central element, or another element
449  // that may or may not be part of the subdomain.
450  const auto& neighbor_mortar_meshes =
451  all_neighbor_mortar_meshes.at(overlap_id);
452  for (const auto& neighbor_mortar_id_and_data :
453  neighbors_mortar_data_.at(overlap_id)) {
454  // No structured bindings because capturing these in lambdas doesn't
455  // work until C++20
456  const auto& neighbor_mortar_id = neighbor_mortar_id_and_data.first;
457  const auto& neighbor_mortar_data = neighbor_mortar_id_and_data.second;
458  const auto& neighbor_direction = neighbor_mortar_id.first;
459  const auto& neighbors_neighbor_id = neighbor_mortar_id.second;
460  // No need to do anything on external boundaries
461  if (neighbors_neighbor_id == ElementId<Dim>::external_boundary_id()) {
462  continue;
463  }
464  const auto& neighbor_orientation =
465  neighbor.neighbors().at(neighbor_direction).orientation();
466  const auto neighbors_neighbor_direction =
467  neighbor_orientation(neighbor_direction.opposite());
468  const ::dg::MortarId<Dim> mortar_id_from_neighbors_neighbor{
469  neighbors_neighbor_direction, neighbor_id};
470  const auto send_mortar_data =
471  [&neighbor_orientation, &neighbor_mortar_meshes,
472  &neighbor_mortar_data, &neighbor_mortar_id,
473  &neighbor_direction](auto& remote_mortar_data) noexcept {
474  const auto& neighbor_mortar_mesh =
475  neighbor_mortar_meshes.at(neighbor_mortar_id);
476  auto oriented_neighbor_mortar_data =
477  neighbor_mortar_data.local_data(temporal_id);
478  if (not neighbor_orientation.is_aligned()) {
479  oriented_neighbor_mortar_data.orient_on_slice(
480  neighbor_mortar_mesh.extents(),
481  neighbor_direction.dimension(), neighbor_orientation);
482  }
483  remote_mortar_data.remote_insert(
484  temporal_id, std::move(oriented_neighbor_mortar_data));
485  };
486  if (neighbors_neighbor_id == central_element.id() and
487  mortar_id_from_neighbors_neighbor == mortar_id) {
488  send_mortar_data(central_mortar_data_.at(mortar_id));
489  continue;
490  }
491  // Determine whether the neighbor's neighbor overlaps with the
492  // subdomain and find its overlap ID if it does.
493  const auto neighbors_neighbor_overlap_id =
494  [&local_all_neighbor_mortar_meshes = all_neighbor_mortar_meshes,
495  &neighbors_neighbor_id,
496  &mortar_id_from_neighbors_neighbor]() noexcept
498  for (const auto& [local_overlap_id, local_mortar_meshes] :
499  local_all_neighbor_mortar_meshes) {
500  if (local_overlap_id.second != neighbors_neighbor_id) {
501  continue;
502  }
503  for (const auto& local_mortar_id_and_mesh : local_mortar_meshes) {
504  if (local_mortar_id_and_mesh.first ==
505  mortar_id_from_neighbors_neighbor) {
506  return local_overlap_id;
507  }
508  }
509  }
510  return std::nullopt;
511  }();
512  if (neighbors_neighbor_overlap_id.has_value()) {
513  // The neighbor's neighbor is part of the subdomain so we copy the
514  // mortar data over. Once the loop is complete we will also have
515  // received mortar data back. At that point, both neighbors have a
516  // copy of each other's mortar data, which is the subject of the
517  // possible optimizations mentioned above. Note that the data may
518  // differ by orientations.
519  send_mortar_data(
520  neighbors_mortar_data_[*neighbors_neighbor_overlap_id]
521  [mortar_id_from_neighbors_neighbor]);
522  } else {
523  // The neighbor's neighbor does not overlap with the subdomain, so
524  // we don't copy mortar data and also don't expect to receive any.
525  // Instead, we assume the data on it is zero and manufacture
526  // appropriate remote boundary data.
527  const auto& neighbors_neighbor_mortar_mesh =
528  all_neighbors_neighbor_mortar_meshes.at(overlap_id)
529  .at(neighbor_mortar_id);
530  auto zero_mortar_data = elliptic::dg::zero_boundary_data_on_mortar<
531  typename System::primal_fields, typename System::primal_fluxes>(
532  neighbors_neighbor_direction,
533  all_neighbors_neighbor_meshes.at(overlap_id)
534  .at(neighbor_mortar_id),
535  all_neighbors_neighbor_face_normal_magnitudes.at(overlap_id)
536  .at(neighbor_mortar_id),
537  neighbors_neighbor_mortar_mesh,
538  all_neighbors_neighbor_mortar_sizes.at(overlap_id)
539  .at(neighbor_mortar_id));
540  // The data is zero, but auxiliary quantities such as the face
541  // normal magnitude may need re-orientation
542  if (not neighbor_orientation.is_aligned()) {
543  zero_mortar_data.orient_on_slice(
544  neighbors_neighbor_mortar_mesh.extents(),
545  neighbors_neighbor_direction.dimension(),
546  neighbor_orientation.inverse_map());
547  }
548  neighbors_mortar_data_.at(overlap_id)
549  .at(neighbor_mortar_id)
550  .remote_insert(temporal_id, std::move(zero_mortar_data));
551  }
552  } // loop over neighbor's mortars
553  } // loop over neighbors
554  } // loop over directions
555 
556  // 2. Apply the operator on all elements in the subdomain
557  //
558  // Apply on central element
559  db::apply<apply_args_tags>(
560  [this, &result, &operand](const auto&... args) noexcept {
561  elliptic::dg::apply_operator<System, linearized>(
562  make_not_null(&result->element_data),
563  make_not_null(&central_mortar_data_), operand.element_data,
564  central_primal_fluxes_, args...);
565  },
566  box, temporal_id, sources_args);
567  // Apply on neighbors
568  for (const auto& [direction, neighbors] : central_element.neighbors()) {
569  const auto& orientation = neighbors.orientation();
570  const auto direction_from_neighbor = orientation(direction.opposite());
571  for (const auto& neighbor_id : neighbors) {
572  const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
573  neighbor_id};
574  const auto& overlap_extent = all_overlap_extents.at(overlap_id);
575  const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
576 
577  if (UNLIKELY(overlap_extent == 0)) {
578  continue;
579  }
580 
581  elliptic::util::apply_at<apply_args_tags_overlap,
582  args_tags_from_center>(
583  [this, &overlap_id](const auto&... args) noexcept {
584  elliptic::dg::apply_operator<System, linearized>(
585  make_not_null(&extended_results_[overlap_id]),
586  make_not_null(&neighbors_mortar_data_.at(overlap_id)),
587  extended_operand_vars_.at(overlap_id),
588  neighbors_primal_fluxes_.at(overlap_id), args...);
589  },
590  box, overlap_id, temporal_id,
591  elliptic::util::apply_at<sources_args_tags_overlap,
592  args_tags_from_center>(get_items, box,
593  overlap_id));
594 
595  // Restrict the extended operator data back to the subdomain, assuming
596  // we can discard any data outside the overlaps. WARNING: This
597  // assumption may break with changes to the DG operator that affect its
598  // sparsity. For example, multiplying the DG operator with the _full_
599  // inverse mass-matrix ("massless" scheme with no "mass-lumping"
600  // approximation) means that lifted boundary corrections bleed into the
601  // volume.
602  if (UNLIKELY(
603  result->overlap_data[overlap_id].number_of_grid_points() !=
604  operand.overlap_data.at(overlap_id).number_of_grid_points())) {
605  result->overlap_data[overlap_id].initialize(
606  operand.overlap_data.at(overlap_id).number_of_grid_points());
607  }
609  make_not_null(&result->overlap_data[overlap_id]),
610  extended_results_.at(overlap_id), neighbor_mesh.extents(),
611  overlap_extent, direction_from_neighbor);
612  } // loop over neighbors
613  } // loop over directions
614  }
615 
616  private:
617  Variables<typename System::auxiliary_fields> central_auxiliary_vars_{};
618  Variables<typename System::primal_fluxes> central_primal_fluxes_{};
619  Variables<typename System::auxiliary_fluxes> central_auxiliary_fluxes_{};
621  Dim, Variables<typename System::auxiliary_fields>>
622  neighbors_auxiliary_vars_{};
624  Variables<typename System::primal_fluxes>>
625  neighbors_primal_fluxes_{};
627  Dim, Variables<typename System::auxiliary_fluxes>>
628  neighbors_auxiliary_fluxes_{};
630  Variables<typename System::primal_fields>>
631  extended_operand_vars_{};
633  Variables<typename System::primal_fields>>
634  extended_results_{};
636  Dim, elliptic::dg::MortarData<size_t, typename System::primal_fields,
637  typename System::primal_fluxes>>
638  central_mortar_data_{};
641  size_t, typename System::primal_fields,
642  typename System::primal_fluxes>>>
643  neighbors_mortar_data_{};
644 };
645 
646 } // namespace elliptic::dg::subdomain_operator
domain::Tags::UnnormalizedFaceNormal
Definition: FaceNormal.hpp:112
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
std::false_type
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
std::pair
Tags::MortarSize
Definition: Tags.hpp:48
FaceNormal.hpp
Tags.hpp
domain::Tags::Element
Definition: Tags.hpp:97
LinearSolver::Schwarz::SubdomainOperator
Abstract base class for the subdomain operator, i.e. the linear operator restricted to an element-cen...
Definition: SubdomainOperator.hpp:50
Domain.hpp
LinearSolver::Schwarz::ElementCenteredSubdomainData
Data on an element-centered subdomain.
Definition: ElementCenteredSubdomainData.hpp:35
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
tuple
FixedHashMap::emplace
std::pair< iterator, bool > emplace(Args &&... args) noexcept
Inserts the element if it does not exists.
Definition: FixedHashMap.hpp:141
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:298
Direction< Dim >
elliptic::dg::subdomain_operator::Tags::ExtrudingExtent
The number of points an element-centered subdomain extends into the neighbor, i.e....
Definition: Tags.hpp:21
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
ElementId< Dim >
ElementId.hpp
elliptic::dg::MortarData
Boundary data on both sides of a mortar.
Definition: DgOperator.hpp:215
make_interface_tag
Definition: InterfaceHelpers.hpp:35
Tags::Normalized
Definition: Magnitude.hpp:137
DataBox.hpp
cstddef
Assert.hpp
domain::Tags::InverseJacobian< Dim, Frame::Logical, Frame::Inertial >
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
LinearSolver::Schwarz::data_on_overlap
void data_on_overlap(const gsl::not_null< Tensor< DataType, TensorStructure... > * > restricted_tensor, const Tensor< DataType, TensorStructure... > &tensor, const Index< Dim > &volume_extents, const size_t overlap_extent, const Direction< Dim > &direction) noexcept
The part of the tensor data that lies within the overlap region.
Definition: OverlapHelpers.hpp:139
elliptic::dg::Tags::PenaltyParameter
The prefactor to the penalty term of the numerical flux.
Definition: Tags.hpp:54
elliptic::dg::subdomain_operator
Items related to the restriction of the DG operator to an element-centered subdomain.
Definition: InitializeSubdomain.hpp:60
LinearSolver::Schwarz::extended_overlap_data
void extended_overlap_data(const gsl::not_null< Variables< ExtendedTagsList > * > extended_data, const Variables< OverlapTagsList > &overlap_data, const Index< Dim > &volume_extents, const size_t overlap_extent, const Direction< Dim > &direction) noexcept
Extend the overlap data to the full mesh by filling it with zeros outside the overlap region.
Definition: OverlapHelpers.hpp:253
std::decay_t
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
LinearSolver::Schwarz::Tags::Overlaps
The Tag on the overlap region with each neighbor, i.e. on a region extruding from the central element...
Definition: Tags.hpp:115
Variables.hpp
DirectionMap
Definition: DirectionMap.hpp:15
ApplyAt.hpp
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:367
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
domain::Tags::InternalDirections
Definition: Tags.hpp:270
elliptic::dg::zero_boundary_data_on_mortar
BoundaryData< PrimalMortarFields, PrimalMortarFluxes > zero_boundary_data_on_mortar(const Direction< Dim > &direction, const Mesh< Dim > &mesh, const Scalar< DataVector > &face_normal_magnitude, const Mesh< Dim - 1 > &mortar_mesh, const ::dg::MortarSize< Dim - 1 > &mortar_size) noexcept
Construct elliptic::dg::BoundaryData assuming the variable data on the element is zero,...
Definition: DgOperator.hpp:190
optional
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
domain::Tags::Domain
Definition: Tags.hpp:55
Direction.hpp
FixedHashMap
A hash table with a compile-time specified maximum size and ability to efficiently handle perfect has...
Definition: FixedHashMap.hpp:82
ostream
elliptic::dg::subdomain_operator::SubdomainOperator
The elliptic DG operator on an element-centered subdomain.
Definition: SubdomainOperator.hpp:98
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:235
std::unordered_map
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:42
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
Tags::Mortars
Definition: Tags.hpp:37