NormalCovectorAndMagnitude.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <memory>
8 #include <optional>
9 
10 #include "DataStructures/DataBox/Tag.hpp"
11 #include "DataStructures/DataVector.hpp"
13 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
18 #include "Domain/Structure/DirectionMap.hpp"
19 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/TMPL.hpp"
22 #include "Utilities/TypeTraits/CreateHasTypeAlias.hpp"
23 
24 namespace evolution::dg::Actions::detail {
25 CREATE_HAS_TYPE_ALIAS(inverse_spatial_metric_tag)
26 CREATE_HAS_TYPE_ALIAS_V(inverse_spatial_metric_tag)
27 
28 struct OneOverNormalVectorMagnitude : db::SimpleTag {
29  using type = Scalar<DataVector>;
30 };
31 
32 template <size_t Dim>
33 struct NormalVector : db::SimpleTag {
34  using type = tnsr::I<DataVector, Dim, Frame::Inertial>;
35 };
36 
37 /*!
38  * \brief Computes the normal covector, magnitude of the unnormalized normal
39  * covector, and in curved spacetimes the normal vector.
40  *
41  * The `fields_on_face` argument is used as a temporary buffer via the
42  * `OneOverNormalVectorMagnitude` tag and also used to return the normalized
43  * normal vector via the `NormalVector<Dim>` tag. The `fields_on_face` argument
44  * also serves as the input for the inverse spatial metric when the spacetime is
45  * curved. The `fields_on_face` argument must have tags `NormalVector<Dim>` and
46  * `OneOverNormalVectorMagnitude`. If the system specifies
47  * `System::inverse_spatial_metric_tag` then this tag must also be in the
48  * Variables.
49  */
50 template <typename System, size_t Dim, typename FieldsOnFaceTags>
51 void unit_normal_vector_and_covector_and_magnitude_impl(
52  const gsl::not_null<Scalar<DataVector>*> face_normal_magnitude,
53  const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
54  unit_normal_covector,
55  const gsl::not_null<Variables<FieldsOnFaceTags>*> fields_on_face,
56  const tnsr::i<DataVector, Dim, Frame::Inertial>&
57  unnormalized_normal_covector) noexcept {
58  if constexpr (has_inverse_spatial_metric_tag_v<System>) {
59  using inverse_spatial_metric_tag =
60  typename System::inverse_spatial_metric_tag;
61  auto& normal_vector = get<NormalVector<Dim>>(*fields_on_face);
62  const auto& inverse_spatial_metric =
63  get<inverse_spatial_metric_tag>(*fields_on_face);
64  // Compute unnormalized normal vector
65  for (size_t i = 0; i < Dim; ++i) {
66  normal_vector.get(i) = inverse_spatial_metric.get(i, 0) *
67  get<0>(unnormalized_normal_covector);
68  for (size_t j = 1; j < Dim; ++j) {
69  normal_vector.get(i) += inverse_spatial_metric.get(i, j) *
70  unnormalized_normal_covector.get(j);
71  }
72  }
73  // Perform normalization
74  dot_product(face_normal_magnitude, normal_vector,
75  unnormalized_normal_covector);
76  get(*face_normal_magnitude) = sqrt(get(*face_normal_magnitude));
77  auto& one_over_normal_vector_magnitude =
78  get<OneOverNormalVectorMagnitude>(*fields_on_face);
79  get(one_over_normal_vector_magnitude) = 1.0 / get(*face_normal_magnitude);
80  for (size_t i = 0; i < Dim; ++i) {
81  unit_normal_covector->get(i) = unnormalized_normal_covector.get(i) *
82  get(one_over_normal_vector_magnitude);
83  normal_vector.get(i) *= get(one_over_normal_vector_magnitude);
84  }
85  } else {
86  magnitude(face_normal_magnitude, unnormalized_normal_covector);
87  auto& one_over_normal_vector_magnitude =
88  get<OneOverNormalVectorMagnitude>(*fields_on_face);
89  get(one_over_normal_vector_magnitude) = 1.0 / get(*face_normal_magnitude);
90  for (size_t i = 0; i < Dim; ++i) {
91  unit_normal_covector->get(i) = unnormalized_normal_covector.get(i) *
92  get(one_over_normal_vector_magnitude);
93  }
94  }
95 }
96 
97 /*!
98  * \brief Computes the normal vector, covector, and their magnitude on the face
99  * in the given direction.
100  *
101  * The normal covector and its magnitude are returned using the
102  * `normal_covector_quantities` function argument. The `fields_on_face` argument
103  * is used as a temporary buffer via the `OneOverNormalVectorMagnitude` tag and
104  * also used to return the normalized normal vector via the `NormalVector<Dim>`
105  * tag. The `fields_on_face` argument also serves as the input for the inverse
106  * spatial metric when the spacetime is curved. The `fields_on_face` argument
107  * must have tags `NormalVector<Dim>` and `OneOverNormalVectorMagnitude`. If the
108  * system specifies `System::inverse_spatial_metric_tag` then this tag must also
109  * be in the Variables.
110  */
111 template <typename System, size_t Dim, typename FieldsOnFaceTags>
112 void unit_normal_vector_and_covector_and_magnitude(
114  Dim, std::optional<Variables<tmpl::list<
117  normal_covector_quantities,
118  const gsl::not_null<Variables<FieldsOnFaceTags>*> fields_on_face,
119  const Direction<Dim>& direction,
120  const std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
121  unnormalized_normal_covectors,
123  moving_mesh_map) noexcept {
124  const bool mesh_is_moving = not moving_mesh_map.is_identity();
125  const auto& unnormalized_normal_covector =
126  unnormalized_normal_covectors.at(direction);
127 
128  if (auto& normal_covector_quantity =
129  normal_covector_quantities->at(direction);
130  detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
131  not normal_covector_quantity.has_value()) {
132  if (not normal_covector_quantity.has_value()) {
133  normal_covector_quantity = Variables<
136  fields_on_face->number_of_grid_points()};
137  }
138  detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
140  &get<evolution::dg::Tags::InternalFace::MagnitudeOfNormal>(
141  *normal_covector_quantity)),
144  *normal_covector_quantity)),
145  fields_on_face, unnormalized_normal_covector);
146  }
147 }
148 } // namespace evolution::dg::Actions::detail
evolution::dg::Tags::InternalFace::MagnitudeOfNormal
The magnitude of the unnormalized normal covector to the interface.
Definition: NormalVectorTags.hpp:19
domain::CoordinateMapBase
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:45
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
evolution::dg::Tags::InternalFace::NormalCovector
The normal covector to the interface.
Definition: NormalVectorTags.hpp:25
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
CoordinateMap.hpp
Direction
Definition: Direction.hpp:23
std::sqrt
T sqrt(T... args)
dot_product
void dot_product(const gsl::not_null< Scalar< DataType > * > dot_product, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean dot product of two vectors or one forms.
Definition: DotProduct.hpp:24
CREATE_HAS_TYPE_ALIAS
#define CREATE_HAS_TYPE_ALIAS(ALIAS_NAME)
Generate a type trait to check if a class has a type alias with a particular name,...
Definition: CreateHasTypeAlias.hpp:27
DotProduct.hpp
cstddef
memory
Variables.hpp
DirectionMap
Definition: DirectionMap.hpp:15
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
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
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
Direction.hpp
std::unordered_map
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13