Line data    Source code

       1           0 : // 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"
      12             : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
      13             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "DataStructures/Variables.hpp"
      16             : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
      17             : #include "Domain/Structure/Direction.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) {
      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<
     112             :                               evolution::dg::Tags::MagnitudeOfNormal,
     113             :                               evolution::dg::Tags::NormalCovector<Dim>>>>>*>
     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,
     119             :     const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>&
     120             :         moving_mesh_map) {
     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,
     132             :                                evolution::dg::Tags::NormalCovector<Dim>>>{
     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)),
     138             :         make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
     139             :             *normal_covector_quantity)),
     140             :         fields_on_face, unnormalized_normal_covector);
     141             :   }
     142             : }
     143             : }  // namespace evolution::dg::Actions::detail