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             : #include <limits>
       9             : #include <memory>
      10             : #include <optional>
      11             : #include <string>
      12             : #include <unordered_map>
      13             : #include <vector>
      14             : 
      15             : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
      16             : #include "Domain/CoordinateMaps/Identity.hpp"
      17             : #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp"
      18             : #include "Domain/Creators/TimeDependence/GenerateCoordinateMap.hpp"
      19             : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
      20             : #include "Domain/Structure/ObjectLabel.hpp"
      21             : #include "Options/Auto.hpp"
      22             : #include "Options/Context.hpp"
      23             : #include "Options/String.hpp"
      24             : #include "Utilities/ErrorHandling/Assert.hpp"
      25             : #include "Utilities/GetOutput.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : 
      28             : /// \cond
      29             : namespace domain::FunctionsOfTime {
      30             : class FunctionOfTime;
      31             : }  // namespace domain::FunctionsOfTime
      32             : namespace domain::CoordinateMaps::TimeDependent {
      33             : class Shape;
      34             : }  // namespace domain::CoordinateMaps::TimeDependent
      35             : namespace domain {
      36             : template <typename SourceFrame, typename TargetFrame, typename... Maps>
      37             : class CoordinateMap;
      38             : }  // namespace domain
      39             : 
      40             : namespace Frame {
      41             : struct Distorted;
      42             : struct Grid;
      43             : struct Inertial;
      44             : }  // namespace Frame
      45             : /// \endcond
      46             : 
      47             : namespace domain::creators::time_dependence {
      48             : /*!
      49             :  * \brief A Shape whose inner surface conforms to a surface of constant
      50             :  * Boyer-Lindquist radius, in Kerr-Schild coordinates as given by
      51             :  * domain::CoordinateMaps::TimeDependent::Shape.
      52             :  *
      53             :  * \details This TimeDependence is suitable for use on a spherical shell,
      54             :  * where LMax is the number of l and m spherical harmonics to use in
      55             :  * approximating the Kerr horizon, of mass `Mass` and spin `Spin`. The value
      56             :  * of the Boyer-Lindquist radius to which the inner surface conforms is given
      57             :  * by the value of `inner_radius`. If the user wants the inner surface of the
      58             :  * Shape to conform to a Kerr horizon for a given mass and spin, `inner_radius`
      59             :  * should be the Boyer-Lindquist radius of the outer horizon.
      60             :  *
      61             :  * The actual shape map that is applied will go from the Grid to the Distorted
      62             :  * frame, and then an identity map will go from the Distorted to Inertial frame.
      63             :  * The reasoning behind this is because in basically all use cases (BBH), the
      64             :  * shape map will go to the Distorted frame only and other maps will go from the
      65             :  * Distorted frame to the Inertial frame.
      66             :  *
      67             :  * \note The quantities stored in the FunctionOfTime are not the
      68             :  * complex spherical-harmonic coefficients \f$\lambda_{lm}(t)\f$, but
      69             :  * instead are the real-valued SPHEREPACK coefficients
      70             :  * \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$ used by Spherepack.  The
      71             :  * relationship between these two sets of coefficients is
      72             :  * \f{align}
      73             :  * a_{l0} & = \sqrt{\frac{2}{\pi}}\lambda_{l0}&\qquad l\geq 0,\\
      74             :  * a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(\lambda_{lm})
      75             :  * &\qquad l\geq 1, m\geq 1, \\
      76             :  * b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(\lambda_{lm})
      77             :  * &\qquad l\geq 1, m\geq 1.
      78             :  * \f}
      79             :  * See domain::CoordinateMaps::TimeDependent::Shape for more details.
      80             :  *
      81             :  * \note To use this time dependence with the `control_system::system::Shape`
      82             :  * control system, you must choose the same \tparam Label that the control
      83             :  * system is using.
      84             :  */
      85             : template <domain::ObjectLabel Label>
      86           1 : class Shape final : public TimeDependence<3> {
      87             :  private:
      88           0 :   using ShapeMap = domain::CoordinateMaps::TimeDependent::Shape;
      89           0 :   using Identity = domain::CoordinateMaps::Identity<3>;
      90           0 :   using GridToInertialMap =
      91             :       detail::generate_coordinate_map_t<Frame::Grid, Frame::Inertial,
      92             :                                         tmpl::list<ShapeMap>>;
      93           0 :   using GridToDistortedMap =
      94             :       detail::generate_coordinate_map_t<Frame::Grid, Frame::Distorted,
      95             :                                         tmpl::list<ShapeMap>>;
      96           0 :   using DistortedToInertialMap =
      97             :       detail::generate_coordinate_map_t<Frame::Distorted, Frame::Inertial,
      98             :                                         tmpl::list<Identity>>;
      99             : 
     100             :  public:
     101           0 :   using maps_list = tmpl::list<
     102             :       domain::CoordinateMap<Frame::Grid, Frame::Inertial, ShapeMap>,
     103             :       domain::CoordinateMap<Frame::Grid, Frame::Distorted, ShapeMap>,
     104             :       domain::CoordinateMap<Frame::Distorted, Frame::Inertial, Identity>>;
     105             : 
     106           0 :   static constexpr size_t mesh_dim = 3;
     107             : 
     108           0 :   static std::string name() { return "Shape"s + get_output(Label); }
     109             : 
     110             :   /// \brief The initial time of the function of time.
     111           1 :   struct InitialTime {
     112           0 :     using type = double;
     113           0 :     static constexpr Options::String help = {
     114             :         "The initial time of the function of time"};
     115             :   };
     116             :   /// \brief The max angular resolution `l` of the Shape.
     117           1 :   struct LMax {
     118           0 :     using type = size_t;
     119           0 :     static constexpr Options::String help = {
     120             :         "The max l value of the Ylms used by the Shape map."};
     121             :   };
     122             :   /// \brief The mass of the Kerr black hole.
     123           1 :   struct Mass {
     124           0 :     using type = double;
     125           0 :     static constexpr Options::String help = {"The mass of the Kerr BH."};
     126             :   };
     127             :   /// \brief The dimensionless spin of the Kerr black hole.
     128           1 :   struct Spin {
     129           0 :     using type = std::array<double, 3>;
     130           0 :     static constexpr Options::String help = {
     131             :         "The dim'less spin of the Kerr BH."};
     132             :   };
     133             :   /// \brief Center for the Shape map
     134           1 :   struct Center {
     135           0 :     using type = std::array<double, 3>;
     136           0 :     static constexpr Options::String help = {"Center for the Shape map."};
     137             :   };
     138             :   /// \brief The inner radius of the Shape map, the radius at which
     139             :   /// to begin applying the map.
     140           1 :   struct InnerRadius {
     141           0 :     using type = double;
     142           0 :     static constexpr Options::String help = {
     143             :         "The inner radius of the Shape map."};
     144             :   };
     145             :   /// \brief The outer radius of the Shape map, beyond which
     146             :   /// it is no longer applied.
     147           1 :   struct OuterRadius {
     148           0 :     using type = double;
     149           0 :     static constexpr Options::String help = {
     150             :         "The outer radius of the Shape map."};
     151             :   };
     152             : 
     153           0 :   using options = tmpl::list<InitialTime, LMax, Mass, Spin, Center, InnerRadius,
     154             :                              OuterRadius>;
     155             : 
     156           0 :   static constexpr Options::String help = {
     157             :       "Creates a Shape that conforms to a Kerr horizon of given mass and "
     158             :       "spin."};
     159             : 
     160           0 :   Shape() = default;
     161           0 :   ~Shape() override = default;
     162           0 :   Shape(const Shape&) = delete;
     163           0 :   Shape(Shape&&) = default;
     164           0 :   Shape& operator=(const Shape&) = delete;
     165           0 :   Shape& operator=(Shape&&) = default;
     166             : 
     167           0 :   Shape(double initial_time, size_t l_max, double mass,
     168             :         std::array<double, 3> spin, std::array<double, 3> center,
     169             :         double inner_radius, double outer_radius,
     170             :         const Options::Context& context = {});
     171             : 
     172           1 :   auto get_clone() const -> std::unique_ptr<TimeDependence<mesh_dim>> override;
     173             : 
     174           1 :   auto block_maps_grid_to_inertial(size_t number_of_blocks) const
     175             :       -> std::vector<std::unique_ptr<domain::CoordinateMapBase<
     176             :           Frame::Grid, Frame::Inertial, mesh_dim>>> override;
     177             : 
     178           1 :   auto block_maps_grid_to_distorted(size_t number_of_blocks) const
     179             :       -> std::vector<std::unique_ptr<domain::CoordinateMapBase<
     180             :           Frame::Grid, Frame::Distorted, mesh_dim>>> override;
     181             : 
     182           1 :   auto block_maps_distorted_to_inertial(size_t number_of_blocks) const
     183             :       -> std::vector<std::unique_ptr<domain::CoordinateMapBase<
     184             :           Frame::Distorted, Frame::Inertial, mesh_dim>>> override;
     185             : 
     186           1 :   auto functions_of_time(const std::unordered_map<std::string, double>&
     187             :                              initial_expiration_times = {}) const
     188             :       -> std::unordered_map<
     189             :           std::string,
     190             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
     191             : 
     192             :  private:
     193             :   template <domain::ObjectLabel OtherLabel>
     194             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     195           0 :   friend bool operator==(const Shape<OtherLabel>& lhs,
     196             :                          const Shape<OtherLabel>& rhs);
     197             : 
     198           0 :   using TransitionFunction = domain::CoordinateMaps::
     199             :       ShapeMapTransitionFunctions::ShapeMapTransitionFunction;
     200             : 
     201           0 :   GridToInertialMap grid_to_inertial_map() const;
     202           0 :   GridToDistortedMap grid_to_distorted_map() const;
     203           0 :   static DistortedToInertialMap distorted_to_inertial_map();
     204             : 
     205           0 :   double initial_time_{std::numeric_limits<double>::signaling_NaN()};
     206           0 :   size_t l_max_{2};
     207           0 :   double mass_{std::numeric_limits<double>::signaling_NaN()};
     208           0 :   std::array<double, 3> spin_{
     209             :       make_array<3>(std::numeric_limits<double>::signaling_NaN())};
     210           0 :   std::array<double, 3> center_{
     211             :       make_array<3>(std::numeric_limits<double>::signaling_NaN())};
     212           0 :   inline static const std::string function_of_time_name_{"Shape" +
     213             :                                                          get_output(Label)};
     214           0 :   double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
     215           0 :   double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
     216           0 :   std::unique_ptr<TransitionFunction> transition_func_;
     217             : };
     218             : 
     219             : template <domain::ObjectLabel Label>
     220           0 : bool operator!=(const Shape<Label>& lhs, const Shape<Label>& rhs);
     221             : }  // namespace domain::creators::time_dependence