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 <optional>
       9             : 
      10             : #include "DataStructures/Tensor/Tensor.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13             : namespace domain::CoordinateMaps {
      14             : 
      15             : /*!
      16             :  * \brief Distorts cartesian coordinates \f$x^i\f$ such that a coordinate sphere
      17             :  * \f$\delta_{ij}x^ix^j=C^2\f$ is mapped to an ellipsoid of constant
      18             :  * Kerr-Schild radius \f$r=C\f$.
      19             :  *
      20             :  * The Kerr-Schild radius \f$r\f$ is defined as the largest positive
      21             :  * root of
      22             :  *
      23             :  *  \f{equation*}{
      24             :  *  r^4 - r^2 (x^2 - a^2) - (\vec{a}\cdot \vec{x})^2 = 0.
      25             :  *  \f}
      26             :  *
      27             :  * In this equation, \f$\vec{x}\f$ are the coordinates of the (distorted)
      28             :  * surface, and \f$\vec{a}=\vec{S}/M\f$ is the spin-parameter of the black hole.
      29             :  * This is equivalent to the implicit definition of \f$r\f$ in Eq. (47) of
      30             :  * \cite Lovelace2008tw. The black hole is assumed to be at the origin.
      31             :  *
      32             :  * \details Given a spin vector \f$\vec{a}\f$, we define the Kerr-Schild radius:
      33             :  *
      34             :  * \f{equation}{r(\vec{x}) = \sqrt{\frac{x^2 - a^2 + \sqrt{(x^2 -
      35             :  * a^2)^2 + 4 (\vec{x} \cdot \vec{a})^2}}{2}} \f}
      36             :  *
      37             :  * We also define the auxiliary variable:
      38             :  *
      39             :  * \f{equation}{s(\vec{x}) = \frac{x^2(x^2 + a^2)}{x^4 + (\vec{x}
      40             :  * \cdot \vec{a})^2} \f}
      41             :  *
      42             :  * The map is then given by:
      43             :  *
      44             :  * \f{equation}{\vec{x}(\vec{\xi}) = \vec{\xi} \sqrt{s(\vec{\xi})} \f}
      45             :  *
      46             :  * The inverse map is given by:
      47             :  * \f{equation}{\vec{\xi}(\vec{x}) = \vec{x} \frac{r(\vec{x})}{|\vec{x}|}
      48             :  * \f}
      49             :  *
      50             :  * The jacobian is:
      51             :  * \f{equation}{ \frac{\partial \xi^i}{\partial x^j} =
      52             :  * \delta_{ij} \sqrt{s(\vec{\xi})} +
      53             :  * \frac{\xi_j}{2 \sqrt{s(\vec{\xi})}} \frac{4\xi^2(1 -
      54             :  * s(\vec{\xi})) \xi_j + 2a^2\xi_i + 2 s(\vec{\xi})
      55             :  * (\vec{\xi} \cdot \vec{a}) a_i}{\xi^4 + (\vec{\xi} \cdot
      56             :  * \vec{a})^2} \f}
      57             :  *
      58             :  * The inverse jacobian is:
      59             :  * \f{equation}{\frac{\partial x^i}{\partial \xi^j} =
      60             :  * \delta_{ij}\frac{r(\vec{x})}{|\vec{x}|}
      61             :  * + \frac{r(\vec{x})^2 x_i + (\vec{x} \cdot \vec{a})a_i}{2r(\vec{x})^3
      62             :  * - r(\vec{x}) (x^2 - a^2)}\frac{x_j}{|\vec{x}|} - \frac{r(\vec{x})
      63             :  * x_i x_j}{x^3} \f}
      64             :  */
      65           1 : class KerrHorizonConforming {
      66             :  public:
      67           0 :   KerrHorizonConforming() = default;
      68           0 :   static constexpr size_t dim = 3;
      69             :   /*!
      70             :    * \brief Constructs a Kerr horizon conforming map.
      71             :    *
      72             :    * \param mass The Kerr mass parameter $M$
      73             :    * \param dimensionless_spin The dimensionless spin $\vec{\chi} = \vec{a} / M
      74             :    * = \vec{S} / M^2$, where $M$ is the Kerr mass parameter, $\vec{S}$ is the
      75             :    * angular momentum or quasilocal spin, and $\vec{a}$ is the Kerr spin
      76             :    * parameter.
      77             :    *
      78             :    * \note The horizon depends only on the dimensionful spin parameter $\vec{a}
      79             :    * = M \vec{\chi}$. This constructor takes $M$ and $\vec{\chi}$ separately for
      80             :    * consistency with other code such as gr::Solutions::KerrSchild, and hence to
      81             :    * avoid bugs where the wrong spin quantity is used accidentally.
      82             :    */
      83           1 :   explicit KerrHorizonConforming(const double mass,
      84             :                                  std::array<double, 3> dimensionless_spin);
      85             : 
      86             :   template <typename T>
      87           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
      88             : 
      89           0 :   std::optional<std::array<double, 3>> inverse(
      90             :       const std::array<double, 3>& target_coords) const;
      91             : 
      92             :   template <typename T>
      93           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
      94             :       const std::array<T, 3>& source_coords) const;
      95             : 
      96             :   template <typename T>
      97           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
      98             :       const std::array<T, 3>& source_coords) const;
      99             : 
     100           0 :   bool is_identity() const {
     101             :     return spin_parameter_ == std::array<double, 3>{0., 0., 0.};
     102             :   }
     103             : 
     104           0 :   static constexpr bool supports_hessian{false};
     105             : 
     106           0 :   friend bool operator==(const KerrHorizonConforming& lhs,
     107             :                          const KerrHorizonConforming& rhs);
     108             : 
     109           0 :   void pup(PUP::er& p);
     110             : 
     111             :  private:
     112             :   template <typename T>
     113           0 :   void stretch_factor_square(gsl::not_null<T*> result,
     114             :                              const std::array<T, 3>& source_coords) const;
     115             : 
     116           0 :   std::array<double, 3> spin_parameter_;
     117           0 :   double spin_mag_sq_;
     118             : };
     119           0 : bool operator!=(const KerrHorizonConforming& lhs,
     120             :                 const KerrHorizonConforming& rhs);
     121             : 
     122             : }  // namespace domain::CoordinateMaps