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/Actions/ComputeTimeDerivativeHelpers.hpp"
20 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/TMPL.hpp"
23 
24 namespace evolution::dg::Actions::detail {
25 struct OneOverNormalVectorMagnitude : db::SimpleTag {
26  using type = Scalar<DataVector>;
27 };
28 
29 template <size_t Dim>
30 struct NormalVector : db::SimpleTag {
31  using type = tnsr::I<DataVector, Dim, Frame::Inertial>;
32 };
33 
34 /*!
35  * \brief Computes the normal covector, magnitude of the unnormalized normal
36  * covector, and in curved spacetimes the normal vector.
37  *
38  * The `fields_on_face` argument is used as a temporary buffer via the
39  * `OneOverNormalVectorMagnitude` tag and also used to return the normalized
40  * normal vector via the `NormalVector<Dim>` tag. The `fields_on_face` argument
41  * also serves as the input for the inverse spatial metric when the spacetime is
42  * curved. The `fields_on_face` argument must have tags `NormalVector<Dim>` and
43  * `OneOverNormalVectorMagnitude`. If the system specifies
44  * `System::inverse_spatial_metric_tag` then this tag must also be in the
45  * Variables.
46  */
47 template <typename System, size_t Dim, typename FieldsOnFaceTags>
48 void unit_normal_vector_and_covector_and_magnitude_impl(
49  const gsl::not_null<Scalar<DataVector>*> face_normal_magnitude,
50  const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
51  unit_normal_covector,
52  const gsl::not_null<Variables<FieldsOnFaceTags>*> fields_on_face,
53  const tnsr::i<DataVector, Dim, Frame::Inertial>&
54  unnormalized_normal_covector) noexcept {
55  if constexpr (has_inverse_spatial_metric_tag_v<System>) {
56  using inverse_spatial_metric_tag =
57  typename System::inverse_spatial_metric_tag;
58  auto& normal_vector = get<NormalVector<Dim>>(*fields_on_face);
59  const auto& inverse_spatial_metric =
60  get<inverse_spatial_metric_tag>(*fields_on_face);
61  // Compute unnormalized normal vector
62  for (size_t i = 0; i < Dim; ++i) {
63  normal_vector.get(i) = inverse_spatial_metric.get(i, 0) *
64  get<0>(unnormalized_normal_covector);
65  for (size_t j = 1; j < Dim; ++j) {
66  normal_vector.get(i) += inverse_spatial_metric.get(i, j) *
67  unnormalized_normal_covector.get(j);
68  }
69  }
70  // Perform normalization
71  dot_product(face_normal_magnitude, normal_vector,
72  unnormalized_normal_covector);
73  get(*face_normal_magnitude) = sqrt(get(*face_normal_magnitude));
74  auto& one_over_normal_vector_magnitude =
75  get<OneOverNormalVectorMagnitude>(*fields_on_face);
76  get(one_over_normal_vector_magnitude) = 1.0 / get(*face_normal_magnitude);
77  for (size_t i = 0; i < Dim; ++i) {
78  unit_normal_covector->get(i) = unnormalized_normal_covector.get(i) *
79  get(one_over_normal_vector_magnitude);
80  normal_vector.get(i) *= get(one_over_normal_vector_magnitude);
81  }
82  } else {
83  magnitude(face_normal_magnitude, unnormalized_normal_covector);
84  auto& one_over_normal_vector_magnitude =
85  get<OneOverNormalVectorMagnitude>(*fields_on_face);
86  get(one_over_normal_vector_magnitude) = 1.0 / get(*face_normal_magnitude);
87  for (size_t i = 0; i < Dim; ++i) {
88  unit_normal_covector->get(i) = unnormalized_normal_covector.get(i) *
89  get(one_over_normal_vector_magnitude);
90  }
91  }
92 }
93 
94 /*!
95  * \brief Computes the normal vector, covector, and their magnitude on the face
96  * in the given direction.
97  *
98  * The normal covector and its magnitude are returned using the
99  * `normal_covector_quantities` function argument. The `fields_on_face` argument
100  * is used as a temporary buffer via the `OneOverNormalVectorMagnitude` tag and
101  * also used to return the normalized normal vector via the `NormalVector<Dim>`
102  * tag. The `fields_on_face` argument also serves as the input for the inverse
103  * spatial metric when the spacetime is curved. The `fields_on_face` argument
104  * must have tags `NormalVector<Dim>` and `OneOverNormalVectorMagnitude`. If the
105  * system specifies `System::inverse_spatial_metric_tag` then this tag must also
106  * be in the Variables.
107  */
108 template <typename System, size_t Dim, typename FieldsOnFaceTags>
109 void unit_normal_vector_and_covector_and_magnitude(
110  const gsl::not_null<
111  DirectionMap<Dim, std::optional<Variables<tmpl::list<
114  normal_covector_quantities,
115  const gsl::not_null<Variables<FieldsOnFaceTags>*> fields_on_face,
116  const Direction<Dim>& direction,
117  const std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
118  unnormalized_normal_covectors,
120  moving_mesh_map) noexcept {
121  const bool mesh_is_moving = not moving_mesh_map.is_identity();
122  const auto& unnormalized_normal_covector =
123  unnormalized_normal_covectors.at(direction);
124 
125  if (auto& normal_covector_quantity =
126  normal_covector_quantities->at(direction);
127  detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
128  not normal_covector_quantity.has_value()) {
129  if (not normal_covector_quantity.has_value()) {
130  normal_covector_quantity =
131  Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
133  fields_on_face->number_of_grid_points()};
134  }
135  detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
136  make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
137  *normal_covector_quantity)),
139  *normal_covector_quantity)),
140  fields_on_face, unnormalized_normal_covector);
141  }
142 }
143 } // namespace evolution::dg::Actions::detail
domain::CoordinateMapBase
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:45
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:83
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
CoordinateMap.hpp
Direction< Dim >
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:25
DotProduct.hpp
cstddef
evolution::dg::Tags::MagnitudeOfNormal
The magnitude of the unnormalized normal covector to the interface.
Definition: NormalVectorTags.hpp:18
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
evolution::dg::Tags::NormalCovector
The normal covector to the interface.
Definition: NormalVectorTags.hpp:24
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: ReadSpecPiecewisePolynomial.hpp:13