KerrHorizonConforming.hpp
1 // 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 
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\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  */
67  public:
68  KerrHorizonConforming() = default;
69  static constexpr size_t dim = 3;
70  /*!
71  * constructs a Kerr horizon conforming map.
72  * \param spin : the dimensionless spin
73  */
75 
76  template <typename T>
78  const std::array<T, 3>& source_coords) const noexcept;
79 
81  const std::array<double, 3>& target_coords) const noexcept;
82 
83  template <typename T>
84  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
85  const std::array<T, 3>& source_coords) const noexcept;
86 
87  template <typename T>
88  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
89  const std::array<T, 3>& source_coords) const noexcept;
90 
91  bool is_identity() const noexcept {
92  return spin_ == std::array<double, 3>{0., 0., 0.};
93  }
94 
95  friend bool operator==(const KerrHorizonConforming& lhs,
96  const KerrHorizonConforming& rhs) noexcept;
97 
98  void pup(PUP::er& p);
99 
100  private:
101  template <typename T>
102  void stretch_factor_square(
103  const gsl::not_null<tt::remove_cvref_wrap_t<T>*> result,
104  const std::array<T, 3>& source_coords) const noexcept;
105 
106  std::array<double, 3> spin_;
107  double spin_mag_sq_;
108 };
109 bool operator!=(const KerrHorizonConforming& lhs,
110  const KerrHorizonConforming& rhs) noexcept;
111 
112 } // namespace domain::CoordinateMaps
domain::CoordinateMaps::KerrHorizonConforming
Distorts cartesian coordinates such that a coordinate sphere is mapped to an ellipsoid of constant ...
Definition: KerrHorizonConforming.hpp:66
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
cstddef
array
domain::CoordinateMaps
Contains all coordinate maps.
Definition: Affine.hpp:23
Gsl.hpp
Tensor.hpp
optional
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13