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