Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <string>
       7             : 
       8             : #include "DataStructures/DataBox/Tag.hpp"
       9             : #include "DataStructures/DataBox/TagName.hpp"
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : #include "Utilities/SetNumberOfGridPoints.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// @{
      18             : /*!
      19             :  * \ingroup TensorGroup
      20             :  * \brief Compute the Euclidean magnitude of a rank-1 tensor
      21             :  *
      22             :  * \details
      23             :  * Computes the square root of the sum of the squares of the components of
      24             :  * the rank-1 tensor.
      25             :  */
      26             : template <typename DataType, typename Index>
      27           1 : Scalar<DataType> magnitude(
      28             :     const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector) {
      29             :   return Scalar<DataType>{sqrt(get(dot_product(vector, vector)))};
      30             : }
      31             : 
      32             : template <typename DataType, typename Index>
      33           1 : void magnitude(const gsl::not_null<Scalar<DataType>*> magnitude,
      34             :                const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector) {
      35             :   set_number_of_grid_points(magnitude, vector);
      36             :   dot_product(magnitude, vector, vector);
      37             :   get(*magnitude) = sqrt(get(*magnitude));
      38             : }
      39             : /// @}
      40             : 
      41             : /// @{
      42             : /*!
      43             :  * \ingroup TensorGroup
      44             :  * \brief Compute the magnitude of a rank-1 tensor
      45             :  *
      46             :  * \details
      47             :  * Returns the square root of the input tensor contracted twice with the given
      48             :  * metric. The absolute value is taken of the dot product to ensure a real
      49             :  * result for spacetime (co-)vectors.
      50             :  */
      51             : template <typename DataType, typename Index>
      52           1 : Scalar<DataType> magnitude(
      53             :     const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector,
      54             :     const Tensor<DataType, Symmetry<1, 1>,
      55             :                  index_list<change_index_up_lo<Index>,
      56             :                             change_index_up_lo<Index>>>& metric) {
      57             :   Scalar<DataType> local_magnitude{get_size(get<0>(vector))};
      58             :   magnitude(make_not_null(&local_magnitude), vector, metric);
      59             :   return local_magnitude;
      60             : }
      61             : 
      62             : template <typename DataType, typename Index>
      63           1 : void magnitude(const gsl::not_null<Scalar<DataType>*> magnitude,
      64             :                const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector,
      65             :                const Tensor<DataType, Symmetry<1, 1>,
      66             :                             index_list<change_index_up_lo<Index>,
      67             :                                        change_index_up_lo<Index>>>& metric) {
      68             :   dot_product(magnitude, vector, vector, metric);
      69             :   if constexpr (Index::index_type == IndexType::Spacetime) {
      70             :     // Spacetime dot product can be negative, so take the absolute value to
      71             :     // define the magnitude
      72             :     get(*magnitude) = sqrt(abs(get(*magnitude)));
      73             :   } else {
      74             :     get(*magnitude) = sqrt(get(*magnitude));
      75             :   }
      76             : }
      77             : /// @}
      78             : 
      79             : namespace Tags {
      80             : /// \ingroup DataBoxTagsGroup
      81             : /// \ingroup DataStructuresGroup
      82             : /// The magnitude of a (co)vector
      83             : template <typename Tag>
      84           1 : struct Magnitude : db::PrefixTag, db::SimpleTag {
      85           0 :   using tag = Tag;
      86           0 :   using type = Scalar<DataVector>;
      87             : };
      88             : 
      89             : /// \ingroup DataBoxTagsGroup
      90             : /// \ingroup DataStructuresGroup
      91             : /// The Euclidean magnitude of a (co)vector
      92             : ///
      93             : /// This tag inherits from `Tags::Magnitude<Tag>`
      94             : template <typename Tag>
      95           1 : struct EuclideanMagnitude : Magnitude<Tag>, db::ComputeTag {
      96           0 :   using base = Magnitude<Tag>;
      97           0 :   using return_type = typename base::type;
      98           0 :   static constexpr auto function =
      99             :       static_cast<void (*)(const gsl::not_null<return_type*>,
     100             :                            const typename Tag::type&)>(&magnitude);
     101           0 :   using argument_tags = tmpl::list<Tag>;
     102             : };
     103             : 
     104             : /// \ingroup DataBoxTagsGroup
     105             : /// \ingroup DataStructuresGroup
     106             : /// The magnitude of a (co)vector with respect to a specific metric
     107             : ///
     108             : /// This tag inherits from `Tags::Magnitude<Tag>`
     109             : template <typename Tag, typename MetricTag>
     110           1 : struct NonEuclideanMagnitude : Magnitude<Tag>, db::ComputeTag {
     111           0 :   using base = Magnitude<Tag>;
     112           0 :   using return_type = typename base::type;
     113           0 :   static constexpr auto function = static_cast<void (*)(
     114             :       const gsl::not_null<return_type*>, const typename Tag::type&,
     115             :       const typename MetricTag::type&)>(&magnitude);
     116           0 :   using argument_tags = tmpl::list<Tag, MetricTag>;
     117             : };
     118             : 
     119             : /// \ingroup DataBoxTagsGroup
     120             : /// \ingroup DataStructuresGroup
     121             : /// The normalized (co)vector represented by Tag
     122             : template <typename Tag>
     123           1 : struct Normalized : db::PrefixTag, db::SimpleTag {
     124           0 :   using tag = Tag;
     125           0 :   using type = typename Tag::type;
     126             : };
     127             : 
     128             : /// \ingroup DataBoxTagsGroup
     129             : /// \ingroup DataStructuresGroup
     130             : /// Normalizes the (co)vector represented by Tag
     131             : ///
     132             : /// This tag inherits from `Tags::Normalized<Tag>`
     133             : template <typename Tag>
     134           1 : struct NormalizedCompute : Normalized<Tag>, db::ComputeTag {
     135           0 :   using base = Normalized<Tag>;
     136           0 :   using return_type = typename base::type;
     137           0 :   static void function(const gsl::not_null<return_type*> normalized_vector,
     138             :                        const typename Tag::type& vector_in,
     139             :                        const typename Magnitude<Tag>::type& magnitude) {
     140             :     *normalized_vector = vector_in;
     141             :     for (size_t d = 0; d < normalized_vector->index_dim(0); ++d) {
     142             :       normalized_vector->get(d) /= get(magnitude);
     143             :     }
     144             :   }
     145           0 :   using argument_tags = tmpl::list<Tag, Magnitude<Tag>>;
     146             : };
     147             : 
     148             : }  // namespace Tags