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