Initialization.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <unordered_map>
9 #include <utility>
10 #include <vector>
11 
12 #include "DataStructures/SliceVariables.hpp"
13 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
16 #include "DataStructures/VariablesTag.hpp"
17 #include "Domain/Domain.hpp"
18 #include "Domain/ElementMap.hpp"
19 #include "Domain/FaceNormal.hpp"
22 #include "Domain/Structure/IndexToSliceAt.hpp"
23 #include "Domain/Tags.hpp"
24 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
25 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
26 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
30 #include "Utilities/Gsl.hpp"
31 #include "Utilities/TMPL.hpp"
32 #include "Utilities/TaggedTuple.hpp"
33 
34 namespace elliptic::dg {
35 
36 /// Initialize the background-independent geometry for the elliptic DG operator
37 template <size_t Dim>
39  using return_tags = tmpl::list<
46  using argument_tags = tmpl::list<domain::Tags::InitialExtents<Dim>,
49  void operator()(
52  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Logical>*> logical_coords,
53  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
55  InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>*>
56  inv_jacobian,
57  gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
58  const std::vector<std::array<size_t, Dim>>& initial_extents,
59  const std::vector<std::array<size_t, Dim>>& initial_refinement,
60  const Domain<Dim>& domain,
61  const ElementId<Dim>& element_id) const noexcept;
62 };
63 
64 /// Initialize the geometry on faces and mortars for the elliptic DG operator
65 ///
66 /// Face normals are not yet normalized, and the face-normal magnitudes are not
67 /// yet computed at all. Invoke `elliptic::dg::NormalizeFaceNormals` for this
68 /// purpose, possibly preceded by `elliptic::dg::InitializeBackground` if the
69 /// system has a background metric.
70 template <size_t Dim>
72  private:
73  // These quantities are currently stored on internal and external faces
74  // separately for compatibility with existing code
75  template <typename DirectionsTag>
76  using face_tags = tmpl::list<
78  domain::Tags::Interface<DirectionsTag,
81  DirectionsTag,
84  DirectionsTag,
86 
87  public:
88  using return_tags = tmpl::flatten<
89  tmpl::list<domain::Tags::InternalDirections<Dim>,
91  face_tags<domain::Tags::InternalDirections<Dim>>,
92  face_tags<domain::Tags::BoundaryDirectionsInterior<Dim>>,
93  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
94  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>>;
95  using argument_tags =
96  tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
98  void operator()(
99  gsl::not_null<std::unordered_set<Direction<Dim>>*> internal_directions,
100  gsl::not_null<std::unordered_set<Direction<Dim>>*> external_directions,
102  face_directions_internal,
104  std::unordered_map<Direction<Dim>, tnsr::I<DataVector, Dim>>*>
105  face_inertial_coords_internal,
107  std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>*>
108  face_normals_internal,
110  face_normal_magnitudes_internal,
112  face_directions_external,
114  std::unordered_map<Direction<Dim>, tnsr::I<DataVector, Dim>>*>
115  face_inertial_coords_external,
117  std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>*>
118  face_normals_external,
120  face_normal_magnitudes_external,
121  gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
123  mortar_sizes,
124  const Mesh<Dim>& mesh, const Element<Dim>& element,
125  const ElementMap<Dim, Frame::Inertial>& element_map,
126  const std::vector<std::array<size_t, Dim>>& initial_extents)
127  const noexcept;
128 };
129 
130 /// Initialize background quantities for the elliptic DG operator, possibly
131 /// including the metric necessary for normalizing face normals
132 template <size_t Dim, typename BackgroundFields>
134  using return_tags = tmpl::list<
140  using argument_tags = tmpl::list<
144 
145  template <typename Background>
146  void operator()(
147  const gsl::not_null<Variables<BackgroundFields>*> background_fields,
148  const gsl::not_null<
149  std::unordered_map<Direction<Dim>, Variables<BackgroundFields>>*>
150  face_background_fields_internal,
151  const gsl::not_null<
152  std::unordered_map<Direction<Dim>, Variables<BackgroundFields>>*>
153  face_background_fields_external,
154  const tnsr::I<DataVector, Dim>& inertial_coords, const Mesh<Dim>& mesh,
155  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
156  inv_jacobian,
157  const Element<Dim>& element,
158  const Background& background) const noexcept {
159  *background_fields = variables_from_tagged_tuple(background.variables(
160  inertial_coords, mesh, inv_jacobian, BackgroundFields{}));
161  ASSERT(mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto,
162  "Only Gauss-Lobatto quadrature is currently implemented for "
163  "slicing background fields to faces.");
164  for (const auto& direction : element.internal_boundaries()) {
165  // Possible optimization: Only the background fields in the
166  // System::fluxes_computer::argument_tags are needed on internal faces.
168  make_not_null(&(*face_background_fields_internal)[direction]),
169  *background_fields, mesh.extents(), direction.dimension(),
170  index_to_slice_at(mesh.extents(), direction));
171  }
172  for (const auto& direction : element.external_boundaries()) {
174  make_not_null(&(*face_background_fields_external)[direction]),
175  *background_fields, mesh.extents(), direction.dimension(),
176  index_to_slice_at(mesh.extents(), direction));
177  }
178  }
179 };
180 
181 /// Normalize face normals, possibly using a background metric. Set
182 /// `InvMetricTag` to `void` to normalize face normals with the Euclidean
183 /// magnitude.
184 template <size_t Dim, typename InvMetricTag>
186  using return_tags =
187  tmpl::list<::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>,
189  using argument_tags =
190  tmpl::conditional_t<std::is_same_v<InvMetricTag, void>, tmpl::list<>,
191  tmpl::list<InvMetricTag>>;
192 
193  template <typename... InvMetric>
194  void operator()(
195  const gsl::not_null<tnsr::i<DataVector, Dim>*> face_normal,
196  const gsl::not_null<Scalar<DataVector>*> face_normal_magnitude,
197  const InvMetric&... inv_metric) const noexcept {
198  magnitude(face_normal_magnitude, *face_normal, inv_metric...);
199  for (size_t d = 0; d < Dim; ++d) {
200  face_normal->get(d) /= get(*face_normal_magnitude);
201  }
202  }
203 };
204 
205 } // namespace elliptic::dg
variables_from_tagged_tuple
Variables< tmpl::list< Tags... > > variables_from_tagged_tuple(const tuples::TaggedTuple< Tags... > &tuple) noexcept
Definition: Variables.hpp:817
utility
elliptic::dg::NormalizeFaceNormal
Normalize face normals, possibly using a background metric. Set InvMetricTag to void to normalize fac...
Definition: Initialization.hpp:185
domain::Tags::Coordinates
Definition: Tags.hpp:130
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::unordered_set
Tags::MortarSize
Definition: Tags.hpp:48
FaceNormal.hpp
Tags.hpp
vector
domain::Tags::Element
Definition: Tags.hpp:97
Domain.hpp
Tags::Variables
Definition: VariablesTag.hpp:21
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
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 >
Element::internal_boundaries
const std::unordered_set< Direction< VolumeDim > > & internal_boundaries() const noexcept
The directions of the faces of the Element that are internal boundaries.
Definition: Element.hpp:57
Element
Definition: Element.hpp:29
ElementId< Dim >
elliptic::dg::InitializeBackground
Initialize background quantities for the elliptic DG operator, possibly including the metric necessar...
Definition: Initialization.hpp:133
Tags::Normalized
Definition: Magnitude.hpp:137
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
Assert.hpp
domain::Tags::InverseJacobian< Dim, Frame::Logical, Frame::Inertial >
array
ElementMap
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:33
Domain
A wrapper around a vector of Blocks that represent the computational domain.
Definition: Domain.hpp:40
elliptic::dg
Functionality related to discontinuous Galerkin discretizations of elliptic equations.
Definition: ApplyOperator.hpp:44
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Element::external_boundaries
const std::unordered_set< Direction< VolumeDim > > & external_boundaries() const noexcept
The directions of the faces of the Element that are external boundaries.
Definition: Element.hpp:51
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
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
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
domain::Tags::InitialRefinementLevels
Definition: Tags.hpp:84
Tensor.hpp
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
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
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:235
unordered_map
elliptic::dg::InitializeGeometry
Initialize the background-independent geometry for the elliptic DG operator.
Definition: Initialization.hpp:38
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
domain::Tags::ElementMap
Definition: Tags.hpp:115