Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : 
       9             : #include "DataStructures/DataBox/DataBox.hpp"
      10             : #include "DataStructures/DataBox/Tag.hpp"
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "DataStructures/Transpose.hpp"
      13             : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
      14             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      15             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      16             : #include "Options/String.hpp"
      17             : #include "Parallel/GlobalCache.hpp"
      18             : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
      19             : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
      20             : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
      21             : #include "ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp"
      22             : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
      23             : #include "Utilities/PrettyType.hpp"
      24             : #include "Utilities/ProtocolHelpers.hpp"
      25             : #include "Utilities/Requires.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : #include "Utilities/TaggedTuple.hpp"
      28             : 
      29             : /// \cond
      30             : class DataVector;
      31             : namespace PUP {
      32             : class er;
      33             : }  // namespace PUP
      34             : namespace db {
      35             : template <typename TagsList>
      36             : class DataBox;
      37             : }  // namespace db
      38             : namespace intrp {
      39             : namespace Tags {
      40             : template <typename TemporalId>
      41             : struct TemporalIds;
      42             : }  // namespace Tags
      43             : }  // namespace intrp
      44             : /// \endcond
      45             : 
      46             : namespace intrp {
      47             : 
      48             : namespace OptionHolders {
      49             : /// A surface that conforms to the horizon of a Kerr black hole in
      50             : /// Kerr-Schild coordinates.
      51             : ///
      52             : /// \details In addition to the parameters for the Kerr black hole, this holder
      53             : /// contains the `LMax` which encodes the angular resolution of the spherical
      54             : /// harmonic basis and `AngularOrdering` which encodes the collocation
      55             : /// ordering.
      56           1 : struct KerrHorizon {
      57           0 :   struct LMax {
      58           0 :     using type = size_t;
      59           0 :     static constexpr Options::String help = {
      60             :         "KerrHorizon is expanded in Ylms up to l=LMax"};
      61             :   };
      62           0 :   struct Center {
      63           0 :     using type = std::array<double, 3>;
      64           0 :     static constexpr Options::String help = {"Center of black hole"};
      65             :   };
      66           0 :   struct Mass {
      67           0 :     using type = double;
      68           0 :     static constexpr Options::String help = {"Mass of black hole"};
      69             :   };
      70           0 :   struct DimensionlessSpin {
      71           0 :     using type = std::array<double, 3>;
      72           0 :     static constexpr Options::String help = {
      73             :         "Dimensionless spin of black hole"};
      74             :   };
      75           0 :   struct AngularOrdering {
      76           0 :     using type = intrp::AngularOrdering;
      77           0 :     static constexpr Options::String help = {
      78             :         "Chooses theta,phi ordering in 2d array"};
      79             :   };
      80           0 :   using options =
      81             :       tmpl::list<LMax, Center, Mass, DimensionlessSpin, AngularOrdering>;
      82           0 :   static constexpr Options::String help = {
      83             :       "A Strahlkorper conforming to the horizon (in Kerr-Schild coordinates)"
      84             :       " of a Kerr black hole with a specified center, mass, and spin."};
      85             : 
      86           0 :   KerrHorizon(size_t l_max_in, std::array<double, 3> center_in, double mass_in,
      87             :               std::array<double, 3> dimensionless_spin_in,
      88             :               intrp::AngularOrdering angular_ordering_in,
      89             :               const Options::Context& context = {});
      90             : 
      91           0 :   KerrHorizon() = default;
      92           0 :   KerrHorizon(const KerrHorizon& /*rhs*/) = default;
      93           0 :   KerrHorizon& operator=(const KerrHorizon& /*rhs*/) = delete;
      94           0 :   KerrHorizon(KerrHorizon&& /*rhs*/) = default;
      95           0 :   KerrHorizon& operator=(KerrHorizon&& /*rhs*/) = default;
      96           0 :   ~KerrHorizon() = default;
      97             : 
      98             :   // NOLINTNEXTLINE(google-runtime-references)
      99           0 :   void pup(PUP::er& p);
     100             : 
     101           0 :   size_t l_max{};
     102           0 :   std::array<double, 3> center{};
     103           0 :   double mass{};
     104           0 :   std::array<double, 3> dimensionless_spin{};
     105           0 :   intrp::AngularOrdering angular_ordering;
     106             : };
     107             : 
     108           0 : bool operator==(const KerrHorizon& lhs, const KerrHorizon& rhs);
     109           0 : bool operator!=(const KerrHorizon& lhs, const KerrHorizon& rhs);
     110             : 
     111             : }  // namespace OptionHolders
     112             : 
     113             : namespace OptionTags {
     114             : template <typename InterpolationTargetTag>
     115           0 : struct KerrHorizon {
     116           0 :   using type = OptionHolders::KerrHorizon;
     117           0 :   static constexpr Options::String help{
     118             :       "Options for interpolation onto Kerr horizon."};
     119           0 :   static std::string name() {
     120             :     return pretty_type::name<InterpolationTargetTag>();
     121             :   }
     122           0 :   using group = InterpolationTargets;
     123             : };
     124             : }  // namespace OptionTags
     125             : 
     126             : namespace Tags {
     127             : template <typename InterpolationTargetTag>
     128           0 : struct KerrHorizon : db::SimpleTag {
     129           0 :   using type = OptionHolders::KerrHorizon;
     130           0 :   using option_tags =
     131             :       tmpl::list<OptionTags::KerrHorizon<InterpolationTargetTag>>;
     132             : 
     133           0 :   static constexpr bool pass_metavariables = false;
     134           0 :   static type create_from_options(const type& option) { return option; }
     135             : };
     136             : }  // namespace Tags
     137             : 
     138             : namespace TargetPoints {
     139             : /// \brief Computes points on a Kerr horizon.
     140             : ///
     141             : /// The points are such that they conform to the horizon of a Kerr
     142             : /// black hole (in Kerr-Schild coordinates) with given center, mass,
     143             : /// and dimensionless spin, as specified in the options.
     144             : ///
     145             : /// \note The returned points are not actually on the horizon;
     146             : /// instead, they are collocation points of a Strahlkorper whose
     147             : /// spectral representation matches the horizon shape up to order
     148             : /// LMax, where LMax is the spherical-harmonic order specified in the
     149             : /// options.  As LMax increases, the returned points converge to the
     150             : /// horizon.
     151             : ///
     152             : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
     153             : ///
     154             : /// For requirements on InterpolationTargetTag, see
     155             : /// intrp::protocols::InterpolationTargetTag
     156             : template <typename InterpolationTargetTag, typename Frame>
     157           1 : struct KerrHorizon : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
     158           0 :   using const_global_cache_tags =
     159             :       tmpl::list<Tags::KerrHorizon<InterpolationTargetTag>>;
     160           0 :   using is_sequential = std::false_type;
     161           0 :   using frame = Frame;
     162             : 
     163           0 :   using simple_tags = typename ylm::Tags::items_tags<Frame>;
     164           0 :   using compute_tags = typename ylm::Tags::compute_items_tags<Frame>;
     165             : 
     166             :   template <typename DbTags, typename Metavariables>
     167           0 :   static void initialize(const gsl::not_null<db::DataBox<DbTags>*> box,
     168             :                          const Parallel::GlobalCache<Metavariables>& cache) {
     169             :     const auto& kerr_horizon =
     170             :         Parallel::get<Tags::KerrHorizon<InterpolationTargetTag>>(cache);
     171             : 
     172             :     // Make a Strahlkorper with the correct shape.
     173             :     ylm::Strahlkorper<Frame> strahlkorper(
     174             :         kerr_horizon.l_max, kerr_horizon.l_max,
     175             :         get(gr::Solutions::kerr_horizon_radius(
     176             :             ::ylm::Spherepack(kerr_horizon.l_max, kerr_horizon.l_max)
     177             :                 .theta_phi_points(),
     178             :             kerr_horizon.mass, kerr_horizon.dimensionless_spin)),
     179             :         kerr_horizon.center);
     180             :     Initialization::mutate_assign<simple_tags>(box, std::move(strahlkorper));
     181             :   }
     182             : 
     183             :   template <typename Metavariables, typename DbTags, typename TemporalId>
     184           0 :   static tnsr::I<DataVector, 3, Frame> points(
     185             :       const db::DataBox<DbTags>& box,
     186             :       const tmpl::type_<Metavariables>& metavariables_v,
     187             :       const TemporalId& /*temporal_id*/) {
     188             :     return points(box, metavariables_v);
     189             :   }
     190             :   template <typename Metavariables, typename DbTags>
     191           0 :   static tnsr::I<DataVector, 3, Frame> points(
     192             :       const db::DataBox<DbTags>& box,
     193             :       const tmpl::type_<Metavariables>& /*meta*/) {
     194             :     const auto& kerr_horizon =
     195             :         db::get<Tags::KerrHorizon<InterpolationTargetTag>>(box);
     196             :     if (kerr_horizon.angular_ordering == intrp::AngularOrdering::Strahlkorper) {
     197             :       return db::get<ylm::Tags::CartesianCoords<Frame>>(box);
     198             :     } else {
     199             :       const auto& strahlkorper = db::get<ylm::Tags::Strahlkorper<Frame>>(box);
     200             :       const auto& coords = db::get<ylm::Tags::CartesianCoords<Frame>>(box);
     201             :       const auto physical_extents =
     202             :           strahlkorper.ylm_spherepack().physical_extents();
     203             :       auto transposed_coords =
     204             :           tnsr::I<DataVector, 3, Frame>(get<0>(coords).size());
     205             :       for (size_t i = 0; i < 3; ++i) {
     206             :         transpose(make_not_null(&transposed_coords.get(i)), coords.get(i),
     207             :                   physical_extents[0], physical_extents[1]);
     208             :       }
     209             :       return transposed_coords;
     210             :     }
     211             :   }
     212             : };
     213             : 
     214             : }  // namespace TargetPoints
     215             : }  // namespace intrp