InitializeSubdomain.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <boost/range/join.hpp>
8 #include <cstddef>
9 #include <tuple>
10 #include <type_traits>
11 #include <unordered_map>
12 #include <utility>
13 #include <vector>
14 
16 #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 #include "DataStructures/DataVector.hpp"
18 #include "DataStructures/SliceVariables.hpp"
19 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
22 #include "DataStructures/VariablesTag.hpp"
23 #include "Domain/ElementMap.hpp"
24 #include "Domain/FaceNormal.hpp"
25 #include "Domain/InterfaceHelpers.hpp"
27 #include "Domain/Structure/CreateInitialMesh.hpp"
29 #include "Domain/Structure/DirectionMap.hpp"
32 #include "Domain/Structure/IndexToSliceAt.hpp"
33 #include "Domain/Tags.hpp"
34 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
35 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
37 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
38 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
41 #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
42 #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
43 #include "Utilities/Gsl.hpp"
44 #include "Utilities/MakeArray.hpp"
45 #include "Utilities/Requires.hpp"
46 #include "Utilities/TMPL.hpp"
47 
48 /// \cond
49 namespace Parallel {
50 template <typename Metavariables>
51 struct GlobalCache;
52 } // namespace Parallel
53 namespace tuples {
54 template <typename... Tags>
55 struct TaggedTuple;
56 } // namespace tuples
57 /// \endcond
58 
59 /// Actions related to the DG subdomain operator
61 
62 namespace detail {
63 // Initialize the geometry of a neighbor into which an overlap extends
64 template <size_t Dim>
65 struct InitializeOverlapGeometry {
66  using return_tags = tmpl::list<
73  domain::Tags::Mesh<Dim - 1>, Dim>,
75  ::Tags::MortarSize<Dim - 1>, Dim>>;
76  using argument_tags =
77  tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>>;
78  void operator()(
79  gsl::not_null<size_t*> extruding_extent,
80  gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim>>*> neighbor_meshes,
82  neighbor_face_normal_magnitudes,
84  neighbor_mortar_meshes,
86  neighbor_mortar_sizes,
87  const Element<Dim>& element, const Mesh<Dim>& mesh,
88  const std::vector<std::array<size_t, Dim>>& initial_extents,
89  const ElementId<Dim>& element_id, const Direction<Dim>& overlap_direction,
90  const size_t max_overlap) const noexcept;
91 };
92 } // namespace detail
93 
94 /*!
95  * \brief Initialize the geometry for the DG subdomain operator
96  *
97  * Initializes tags that define the geometry of overlap regions with neighboring
98  * elements. The data needs to be updated if the geometry of neighboring
99  * elements changes.
100  *
101  * Note that the geometry depends on the system and on the choice of background
102  * through the normalization of face normals, which involves a metric.
103  *
104  * DataBox:
105  * - Uses:
106  * - `BackgroundTag`
107  * - `domain::Tags::Element<Dim>`
108  * - `domain::Tags::InitialExtents<Dim>`
109  * - `domain::Tags::InitialRefinementLevels<Dim>`
110  * - `domain::Tags::Domain<Dim>`
111  * - `LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>`
112  * - Adds: Many tags prefixed with `LinearSolver::Schwarz::Tags::Overlaps`. See
113  * `elliptic::dg::Actions::InitializeDomain` and
114  * `elliptic::dg::Actions::initialize_operator` for a complete list.
115  */
116 template <typename System, typename BackgroundTag, typename OptionsGroup>
118  private:
119  static constexpr size_t Dim = System::volume_dim;
120  static constexpr bool is_curved =
121  not std::is_same_v<typename System::inv_metric_tag, void>;
122  static constexpr bool has_background_fields =
123  not std::is_same_v<typename System::background_fields, tmpl::list<>>;
124 
126  using InitializeOverlapGeometry = detail::InitializeOverlapGeometry<Dim>;
129  using NormalizeFaceNormal =
131 
132  template <typename Tag>
133  using overlaps_tag =
135  // Only slice those background fields to internal boundaries that are
136  // necessary for the DG operator, i.e. the background fields in the
137  // System::fluxes_computer::argument_tags
138  using fluxes_non_background_args =
139  tmpl::list_difference<typename System::fluxes_computer::argument_tags,
140  typename System::background_fields>;
141  using background_fields_internal =
142  tmpl::list_difference<typename System::fluxes_computer::argument_tags,
143  fluxes_non_background_args>;
144  // Slice all background fields to external boundaries for use in boundary
145  // conditions
146  using background_fields_external = typename System::background_fields;
147 
148  public:
149  using initialization_tags =
150  tmpl::list<domain::Tags::InitialExtents<Dim>,
152  using const_global_cache_tags =
153  tmpl::list<LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>>;
154  using simple_tags = db::wrap_tags_in<
155  overlaps_tag,
156  tmpl::append<
157  typename InitializeGeometry::return_tags,
158  typename InitializeFacesAndMortars::return_tags,
159  typename InitializeOverlapGeometry::return_tags,
160  tmpl::conditional_t<
161  has_background_fields,
162  tmpl::list<::Tags::Variables<typename System::background_fields>>,
163  tmpl::list<>>,
164  make_interface_tags<background_fields_internal,
166  make_interface_tags<background_fields_external,
168  using compute_tags = tmpl::list<>;
169 
170  template <typename DataBox, typename... InboxTags, typename Metavariables,
171  typename ActionList, typename ParallelComponent>
172  static std::tuple<DataBox&&> apply(
173  DataBox& box, const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
174  const Parallel::GlobalCache<Metavariables>& /*cache*/,
175  const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
176  const ParallelComponent* const /*meta*/) noexcept {
177  if constexpr (tmpl::all<initialization_tags,
178  tmpl::bind<db::tag_is_retrievable, tmpl::_1,
179  tmpl::pin<DataBox>>>::value) {
180  const auto& element = db::get<domain::Tags::Element<Dim>>(box);
181  for (const auto& [direction, neighbors] : element.neighbors()) {
182  const auto& orientation = neighbors.orientation();
183  const auto direction_from_neighbor = orientation(direction.opposite());
184  for (const auto& neighbor_id : neighbors) {
185  const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
186  neighbor_id};
187  // Initialize background-agnostic geometry on overlaps
188  elliptic::util::mutate_apply_at<
190  typename InitializeGeometry::return_tags>,
191  typename InitializeGeometry::argument_tags,
192  typename InitializeGeometry::argument_tags>(
193  InitializeGeometry{}, make_not_null(&box), overlap_id,
194  neighbor_id);
195  // Initialize faces and mortars on overlaps
196  elliptic::util::mutate_apply_at<
198  typename InitializeFacesAndMortars::return_tags>,
200  overlaps_tag,
201  typename InitializeFacesAndMortars::argument_tags>,
202  tmpl::list<>>(InitializeFacesAndMortars{}, make_not_null(&box),
203  overlap_id,
204  db::get<domain::Tags::InitialExtents<Dim>>(box));
205  // Initialize subdomain-specific tags on overlaps
206  elliptic::util::mutate_apply_at<
208  typename InitializeOverlapGeometry::return_tags>,
210  overlaps_tag,
211  typename InitializeOverlapGeometry::argument_tags>,
212  tmpl::list<>>(
213  InitializeOverlapGeometry{}, make_not_null(&box), overlap_id,
214  db::get<domain::Tags::InitialExtents<Dim>>(box), neighbor_id,
215  direction_from_neighbor,
217  box));
218  // Background fields
219  if constexpr (has_background_fields) {
220  initialize_background_fields(make_not_null(&box), overlap_id);
221  }
222  // Normalize face normals
223  normalize_face_normals(make_not_null(&box), overlap_id);
224  } // neighbors in direction
225  } // directions
226  } else {
227  ERROR(
228  "Dependencies not fulfilled. Did you forget to terminate the phase "
229  "after removing options?");
230  }
231  return {std::move(box)};
232  }
233 
234  private:
235  template <typename DbTagsList>
236  static void initialize_background_fields(
237  const gsl::not_null<db::DataBox<DbTagsList>*> box,
238  const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) noexcept {
239  const auto& background = db::get<BackgroundTag>(*box);
241  face_background_fields{};
242  elliptic::util::mutate_apply_at<
243  tmpl::list<overlaps_tag<
246  overlaps_tag,
247  tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
252  tmpl::list<>>(
253  [&background, &face_background_fields](
255  background_fields,
256  const tnsr::I<DataVector, Dim>& inertial_coords,
257  const Mesh<Dim>& mesh,
258  const InverseJacobian<DataVector, Dim, Frame::Logical,
259  Frame::Inertial>& inv_jacobian,
260  const Element<Dim>& element) noexcept {
261  *background_fields = variables_from_tagged_tuple(
262  background.variables(inertial_coords, mesh, inv_jacobian,
263  typename System::background_fields{}));
264  for (const auto& direction :
265  boost::join(element.internal_boundaries(),
266  element.external_boundaries())) {
267  // Slice the background fields to the face instead of evaluating
268  // them on the face coords to avoid re-computing them, and because
269  // this is also what the DG operator currently does. The result is
270  // the same on Gauss-Lobatto grids, but may need adjusting when
271  // adding support for Gauss grids.
272  data_on_slice(make_not_null(&face_background_fields[direction]),
273  *background_fields, mesh.extents(),
274  direction.dimension(),
275  index_to_slice_at(mesh.extents(), direction));
276  }
277  },
278  box, overlap_id);
279  // Move face background fields into DataBox
280  const auto mutate_assign_face_background_field =
281  [&box, &overlap_id, &face_background_fields](
282  auto tag_v, auto directions_tag_v,
283  const Direction<Dim>& direction) noexcept {
284  using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
285  using directions_tag = std::decay_t<decltype(directions_tag_v)>;
286  db::mutate<
288  box, [&face_background_fields, &overlap_id,
289  &direction](const auto stored_value) noexcept {
290  (*stored_value)[overlap_id][direction] =
291  get<tag>(face_background_fields.at(direction));
292  });
293  };
294  const auto& element =
295  db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
296  for (const auto& direction : element.internal_boundaries()) {
297  tmpl::for_each<background_fields_internal>(
298  [&mutate_assign_face_background_field,
299  &direction](auto tag_v) noexcept {
300  mutate_assign_face_background_field(
301  tag_v, domain::Tags::InternalDirections<Dim>{}, direction);
302  });
303  }
304  for (const auto& direction : element.external_boundaries()) {
305  tmpl::for_each<background_fields_external>(
306  [&mutate_assign_face_background_field,
307  &direction](auto tag_v) noexcept {
308  mutate_assign_face_background_field(
310  direction);
311  });
312  }
313  }
314 
315  template <typename DbTagsList>
316  static void normalize_face_normals(
317  const gsl::not_null<db::DataBox<DbTagsList>*> box,
318  const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) noexcept {
319  // Faces of the overlapped element (internal and external)
320  const auto& element =
321  db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
322  for (const auto& direction : element.internal_boundaries()) {
323  elliptic::util::mutate_apply_at<
325  overlaps_tag,
326  make_interface_tags<typename NormalizeFaceNormal::return_tags,
329  overlaps_tag,
330  make_interface_tags<typename NormalizeFaceNormal::argument_tags,
332  tmpl::list<>>(NormalizeFaceNormal{}, box,
333  std::make_tuple(overlap_id, direction));
334  }
335  for (const auto& direction : element.external_boundaries()) {
336  elliptic::util::mutate_apply_at<
338  make_interface_tags<
339  typename NormalizeFaceNormal::return_tags,
342  make_interface_tags<
343  typename NormalizeFaceNormal::argument_tags,
345  tmpl::list<>>(NormalizeFaceNormal{}, box,
346  std::make_tuple(overlap_id, direction));
347  }
348  // Faces on the other side of the overlapped element's mortars
349  const auto& domain = db::get<domain::Tags::Domain<Dim>>(*box);
350  const auto& neighbor_meshes = db::get<overlaps_tag<
352  domain::Tags::Mesh<Dim>, Dim>>>(*box)
353  .at(overlap_id);
354  for (const auto& [direction, neighbors] : element.neighbors()) {
355  const auto& orientation = neighbors.orientation();
356  const auto direction_from_neighbor = orientation(direction.opposite());
357  for (const auto& neighbor_id : neighbors) {
358  const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
359  const auto neighbor_face_mesh =
360  neighbor_meshes.at(mortar_id).slice_away(
361  direction_from_neighbor.dimension());
362  const auto& neighbor_block = domain.blocks()[neighbor_id.block_id()];
363  ElementMap<Dim, Frame::Inertial> neighbor_element_map{
364  neighbor_id, neighbor_block.stationary_map().get_clone()};
365  const auto neighbor_face_normal = unnormalized_face_normal(
366  neighbor_face_mesh, neighbor_element_map, direction_from_neighbor);
367  using neighbor_face_normal_magnitudes_tag = overlaps_tag<
370  Dim>>;
371  if constexpr (is_curved) {
372  const auto& background = db::get<BackgroundTag>(*box);
373  const auto neighbor_face_inertial_coords =
374  neighbor_element_map(interface_logical_coordinates(
375  neighbor_face_mesh, direction_from_neighbor));
376  const auto inv_metric_on_face =
377  get<typename System::inv_metric_tag>(background.variables(
378  neighbor_face_inertial_coords,
379  tmpl::list<typename System::inv_metric_tag>{}));
380  elliptic::util::mutate_apply_at<
381  tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
382  tmpl::list<>>(
383  [&neighbor_face_normal, &inv_metric_on_face](
384  const auto neighbor_face_normal_magnitude) noexcept {
385  magnitude(neighbor_face_normal_magnitude, neighbor_face_normal,
386  inv_metric_on_face);
387  },
388  box, std::make_tuple(overlap_id, mortar_id));
389  } else {
390  elliptic::util::mutate_apply_at<
391  tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
392  tmpl::list<>>(
393  [&neighbor_face_normal](
394  const auto neighbor_face_normal_magnitude) noexcept {
395  magnitude(neighbor_face_normal_magnitude, neighbor_face_normal);
396  },
397  box, std::make_tuple(overlap_id, mortar_id));
398  }
399  } // neighbors
400  } // internal directions
401  }
402 };
403 
404 } // namespace elliptic::dg::subdomain_operator::Actions
variables_from_tagged_tuple
Variables< tmpl::list< Tags... > > variables_from_tagged_tuple(const tuples::TaggedTuple< Tags... > &tuple) noexcept
Definition: Variables.hpp:817
utility
Frame::Inertial
Definition: IndexType.hpp:44
elliptic::dg::NormalizeFaceNormal
Normalize face normals, possibly using a background metric. Set InvMetricTag to void to normalize fac...
Definition: Initialization.hpp:185
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
index_to_slice_at
size_t index_to_slice_at(const Index< Dim > &extents, const Direction< Dim > &direction, const size_t offset=0) noexcept
Finds the index in the perpendicular dimension of an element boundary.
Definition: IndexToSliceAt.hpp:18
std::pair
Tags::MortarSize
Definition: Tags.hpp:48
FaceNormal.hpp
Tags.hpp
vector
MakeArray.hpp
domain::Tags::Element
Definition: Tags.hpp:97
interface_logical_coordinates
tnsr::I< DataVector, VolumeDim, Frame::Logical > interface_logical_coordinates(const Mesh< VolumeDim - 1 > &mesh, const Direction< VolumeDim > &direction) noexcept
Compute the logical coordinates on a face of an Element.
LinearSolver::Schwarz::Tags::MaxOverlap
Number of points a subdomain can overlap with its neighbor.
Definition: Tags.hpp:63
Tags::Variables
Definition: VariablesTag.hpp:21
elliptic::dg::subdomain_operator::Actions::InitializeSubdomain
Initialize the geometry for the DG subdomain operator.
Definition: InitializeSubdomain.hpp:117
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
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:298
elliptic::dg::InitializeFacesAndMortars
Initialize the geometry on faces and mortars for the elliptic DG operator.
Definition: Initialization.hpp:71
Spectral.hpp
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
Element
Definition: Element.hpp:29
ElementId< Dim >
ElementId.hpp
db::mutate
decltype(auto) mutate(const gsl::not_null< DataBox< TagList > * > box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:641
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
elliptic::dg::subdomain_operator::Tags::NeighborMortars
Data on the neighbor's side of a mortar. Used to store data for elements that do not overlap with the...
Definition: Tags.hpp:29
DataBox.hpp
db::data_on_slice
Variables< tmpl::list< TagsToSlice... > > data_on_slice(const db::DataBox< TagsList > &box, const Index< VolumeDim > &element_extents, const size_t sliced_dim, const size_t fixed_index, tmpl::list< TagsToSlice... >) noexcept
Slices volume Tensors from a DataBox into a Variables
Definition: DataOnSlice.hpp:33
cstddef
domain::Tags::InverseJacobian
The inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:165
db::tag_is_retrievable
tmpl::any< typename DataBoxType::tags_list, std::is_base_of< tmpl::pin< Tag >, tmpl::_1 > > tag_is_retrievable
Definition: DataBox.hpp:58
array
LogicalCoordinates.hpp
ElementMap
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:33
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
std::decay_t
elliptic::dg::subdomain_operator::Actions
Actions related to the DG subdomain operator.
Definition: InitializeSubdomain.hpp:60
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
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
DirectionMap
Definition: DirectionMap.hpp:15
ApplyAt.hpp
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Tags::Magnitude
Definition: Magnitude.hpp:98
Gsl.hpp
domain::Tags::InitialRefinementLevels
Definition: Tags.hpp:84
Tensor.hpp
Requires.hpp
domain::Tags::InternalDirections
Definition: Tags.hpp:270
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
magnitude
Scalar< DataType > magnitude(const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector) noexcept
Compute the Euclidean magnitude of a rank-1 tensor.
Definition: Magnitude.hpp:27
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
unordered_map
elliptic::dg::InitializeGeometry
Initialize the background-independent geometry for the elliptic DG operator.
Definition: Initialization.hpp:38
brigand::list_difference
fold< Sequence2, Sequence1, lazy::remove< _state, _element > > list_difference
Obtain the elements of Sequence1 that are not in Sequence2.
Definition: TMPL.hpp:589
Parallel
Functionality for parallelization.
Definition: ElementReceiveInterpPoints.hpp:13
type_traits
unnormalized_face_normal
void unnormalized_face_normal(gsl::not_null< tnsr::i< DataVector, VolumeDim, TargetFrame > * > result, const Mesh< VolumeDim - 1 > &interface_mesh, const ElementMap< VolumeDim, TargetFrame > &map, const Direction< VolumeDim > &direction) noexcept
Compute the outward grid normal on a face of an Element.
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13